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I. Introduction 



In this paper, I will concentrate on those Fokker-Planck models which are most useful for the study 
of rf-driven currents. 1 I will therefore take the plasma to be azimuthally symmetric about the magnetic 
field and homogeneous (representative of the central portion of a tokamak plasma). The Fokker-Planck 
equation then reduces to an equation in time and two velocity (or momentum) dimensions only. This 
simplified model yields a wealth of interesting physics and furthermore illustrates the main numerical 
problems encountered in more complicated situations. In addition to the collision term, the equation will 
include the effects of externally injected rf power via a quasilinear diffusion term, and a dc electric field. 
(The electric field arises whenever the current is time-varying, e.g., during current ramp-up.) Because 
the wave interacts with very fast electrons, relativistic effects are also considered. In addition, the adjoint 
method for solving for moments of the Fokker-Planck equation is discussed. This method allows for a 
great reduction (by orders of magnitude) in the amount of computer time required. 

The paper is divided into three parts: In the first part of the paper, I give the formulation of the Fokker- 
Planck equation. In sec. II, the Fokker-Planck equation and the coordinate systems are introduced. The 
collision operator and approximations to it are given in sees. Ill and IV. Corresponding expressions for 
the quasilinear diffusion operator are given in sec. V. The next part of the paper describes the numerical 
solution of the equation. Its boundary conditions are considered in sec. VI. Sections VII and VIII describe 
the spatial and temporal differencing of the equation. In sec. IX, we describe techniques for obtaining 
the time asymptotic solution to the equation. The last part of the paper describes the incorporation of 
relativistic effects (sec. X) and the adjoint method for solving the Fokker-Planck equation (sec. XI). 

The numerical methods presented here are those used in the Fokker-Planck code used by the author. 
A word about the lineage of this code is in order: Fokker-Planck codes were developed at Livermore 
by Killeen et ai. 2 ' 3 for the study of mirror-machine plasmas. The latest stage in the development of 
these codes is FPPAC 4 which is a two-dimensional multispecies nonlinear Fokker-Planck package. The 
Livermore code was extensively modified by Winsor and Fallon for the study of runaway electrons, which 
was undertaken by Kulsrud et ai.. 5 This code has been used by the author in various studies of current 
drive beginning with lower hybrid current drive. 6 Over the years several further modifications have been 
made, although the basic structure of the code is the same as that of Kulsrud et aJ. 

The assumption of a homogeneous magnetic field is warranted in the study of rf heating and current 
drive in tokamaks if the rf interacts only with circulating particles. This is often not the case (for exam- 
ple during ion- and electron-cyclotron heating) in which case proper account should be taken of trapped 
particles. This has been done in so-called bounce-averaged codes 7 9 in which the distribution function 
is averaged over the bounce motion of the trapped particles. In tokamaks this leads to a modification of 
the coefficients appearing in the Fokker-Planck equation but the numerical treatment of the equation is 
largely unaltered. In machines with more complicated particle orbits, the distribution may be a multi- 
valued function of the velocity coordinates. This occurs in tandem mirrors where there is more than one 
population of trapped particles. In this case, special techniques are required. 9 

II. Preliminaries 

A. The Fokker-Planck equation 

We write the Fokker-Planck equation for the electrons e as 

^ - £C(/ e ,/,) + V - S w + ■ V/ e = 0, (1) 

s 

where q s and m s are the charge and mass of species s, C(f ai fb) is the collision term for species a 
colliding off species b, the sum extends over all the species of the plasma (typically electrons and ions), 
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S w is the wave (u>)-induced quasilinear flux, and E = Eir\\ is the electric field (assumed to be parallel to 
the magnetic field). The quantity q s carries the sign of the charge, thus q e = —e. The subscripts || and _L 
refer to the directions parallel and perpendicular to the magnetic field. The V = d/dv operator operates 
in velocity space. 

Because collisions in a plasma are primarily due to small-angle scattering, the collision term can be 
written as the divergence of a flux 

C(f a J b ) = -V-S a J\ 



in which case eq. (1) can be expressed as 



where 

is the total flux in velocity space, and 



^ +V ' S = °< (2) 



S — S c + S w + S e 



S C = ]TS^, (3) 

s 

Se = —fe, (4) 



are the collisional (c)- and electric-field (e)-induced electron fluxes. 
From eq. (2) we can derive the conservation laws 



/ f e d 3 v + f S • d 2 A = 0, (5a) 

JV J A 



dt J V • i 

d_ 

di . . 

d f m, f- '" • '-' 

dt 



(5b) 



/ m e v/ e c? 3 v+ / m e vS • d 2 A = / m e Sd 3 v, 
Jv J A Jv 

J ^ff e d\ + J A ^fs- d 2 A = ^ m e v • S d 3 v, (5c) 



where V is some volume in velocity space and A is its boundary. These equations are statements of 
conservation of number, momentum, and energy. 

Typically, two types of terms appear in S: a diffusion term and a friction term 

S = -D- V/e+F/ e . (6) 

The wave term is purely diffusive so that F„ = 0, while the electric field term is nondiffusive: D e = 0, 

F e = (7) 

m e 



B. Coordinate systems 

Because of azimuthal symmetry, f e is independent of <f) the angle about the magnetic field. Two 
coordinate systems suggest themselves: the cylindrical coordinate system (v±,v\\,<f>) and the spherical 
coordinate system (v, 9, <f>); see fig. 1. These are related by 

2 2,2 
V =V ± +Vp 

cos 6 — v\\/v. 
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Both of these coordinate systems are useful. In cylindrical coordinates (assuming azimuthal symmetry) 
eq. (6) gives 



v-s = — -!L v± s ± + -£-s n , 

V± OV± OVu " 



S± 
Sn 



D d/e £i df e 



F±fe, 



Similarly, in spherical coordinates we have 

1 d n 

V-S = ^— v 2 S v 
it ov 



1 d 
v sin 9 d9 



sin 9 So , 



OV V ov 

Sg = -D ev ^ - Dge-^ + F e f e . 
Ov v 09 

Transformations between D and F expressed in the two coordinate systems may be achieved by 



(8a) 
(8b) 
(8c) 

(9a) 
(9b) 
(9c) 



(D ±± \ 




( D vv 


Ami i 


= M 


D v9 


D \\± I 


Dg v 


V Am / 




\D 66 



and 



where 



F, 



sc 
sc 
c 2 
s 
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)■ 






sc 


sc 


-s 2 


c 2 


c 2 - 


-s 2 


—sc - 


-sc 











(10a) 



(10b) 



c2 \ 

— sc 

—sc 

c 2 



M = M 1 
N = N 1 



and here we have abbreviated s = sin 9 and c = cos 9. The collision term is most conveniently expressed 
in spherical coordinates and eqs. (10) allow us to transform this term to cylindrical coordinates. On the 
other hand, the rf and electric field terms are written most naturally in cylindrical coordinates, and this 
equation also enables us to express these terms in spherical coordinates. 

In the case of the collision operator, D and F are given in terms of the gradients of potentials 

D oc VVV>, F oc V<p. 

In cylindrical coordinates (with azimuthal symmetry), the relevant components of D and F are easy to 
calculate — we just take the corresponding derivatives of and <j>. In spherical coordinates we have 



(VVV>)™ = 



d 2 ^ 
dv 2 ' 
1 d 2 ^ 
v dvdO 



1 <9V 



(11a) 
(lib) 
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™» = z-£ + *m> (llc) 

(Vtf)„ = ^, (Hd) 

<™* = \w (lle) 

Several important quantities are given in terms of velocity-space moments of the distribution function. 
Three-dimensional velocity space integrations can be carried out in cylindrical coordinates using 

//>OC />OC 
/(v)d 3 v= / dv± dv^Trv^fiv^vn) (12) 
JO J -oo 

and in spherical coordinates using 

/p OO P 7T 

/(v)d 3 v = dv d9 2irv 2 f(v,9)sm9. (13) 
Jo Jo 

C. Legendre harmonics 

It is sometimes useful to decompose the distribution function and the potentials into Legendre har- 
monics P;(cos 9). We write this as 

OO 

f(v,e) = Y,f {l Hv)Pi(cos6), (14) 

2=0 

where 

/ (0 («) = £ f(v ,9)^(008 9) sin 9 d9. (15) 

The Legendre polynomials may be evaluated on a computer using the recurrence relation 10 

Po(p) = 1, (/+ l)P m (M) - (2/ + 1)mW - m-i(M)- 

D. Definitions 

Finally, we define some of the other quantities that we encounter in this paper. The thermal velocity 
of species s is given by 

v ts = J^-, (16) 
V m s 

where T s is the temperature of species s. The thermal collision frequency for the electrons is 

pe/e 

Vte = T te = —g-, (17) 

v te 

where 

r a/ fc= n b qlq 2 h In A*/ b 



and n s is the number density of species s, e is the dielectric constant of free space, and In A a//fc is the 
Coulomb logarithm. The distributions are normalized so that 



J / s (v) d 3 v = n s . 
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In particular, the Maxwellian distribution is 



3/2 



fsm{v) = n s [ — ) cxp(-±m s v 2 /T s ), 



\2irT s 

^ Us (2^)3/2 ex PH^A4)- (18) 

In discussing the applications to rf current drive, there are two quantities in which we will be inter- 
ested: the electron current density 

J = J q e v\\f e (v)cPv (19) 
and the rf power absorbed per unit volume by the plasma 



P = J m e v • S w tfV. (20) 

The efficiency of rf current drive is usually given as the ratio J/ P. 

We shall use S.I. units throughout this paper except that we will measure temperature in units of 
energy. 

III. Collision Operator 

A. The Landau collision operator 



The collision flux is given by the Landau collision integral 

r,2„2 
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S a /b = _ti£_ lnA a /b f U(u) . ( /.(V) dh^') h<?)dUv) \ ^ 

STre^m a J \ m b dv' m a dv J 



where 



U(u) = ; , U = V-V. 



The formula for the Coulomb logarithm In k a / b is given in text books and the NRL Plasma Formulary. 12 
Because it is so insensitive to plasma parameters, in many cases it is adequate to take it to be a constant 
equal to 15. In any case, it is required that In A a / h — In A b / a . The Landau collision operator conserves 
number, momentum, and energy, i.e., 

J C(f a ,f b )d 3 v = 0, (22a) 
J m a vC{f a , f b ) d 3 v + J m b vC(f b ,f a ) d 3 v - 0, (22b) 
^C(/ , f b ) d 3 v + J ^fc(.f b , f a ) rf 3 v = 0. (22c) 



/ 



There is an error on the order of 1/ In A a / b in the Landau collision operator. However, because it has 
so many "nice" properties — the conservation laws of eqs. (22), an _ff-theorem, etc. — it is customary to 
regard eq. (21) as being exact. 
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B. Rosenbluth potentials 



Equation (21) is the most useful form for the collision operator for analytical work. However, it is not 
in a convenient form for numerical computations. Suppose we represent the distribution functions on an 
N x N grid (assuming azimuthal symmetry). Then evaluation of eq. (21) entails 0(N A ) computations, 
because it entails a two-dimensional integral (over v') which must be carried out at each grid location. 
Fortunately, substantial savings may be realized by using an equivalent representation in terms of Rosen- 
bluth potentials. 1314 Here we use the slightly more convenient notation of Trubnikov. 15 We define two 
potentials 



^(v) = -^ / |v-v'|/ b (v')rfV. 
These are called potentials because they satisfy Poisson's equations in velocity space 

V 2 ^(v) =/&(v), VVft(v) =&(v). 
In terms of these potentials, eq. (21) becomes 

S ^ = -D: /6 V/ a (v) + F«/ 6 /a(v), 

na/6 4^r a /\ 



(23a) 
(23b) 



n b 



-VVV>6(v) 



F a/b = _^Ej_I^V0 6 (v). 



(24a) 
(24b) 

(24c) 



n b m b 

(An equivalent form of this equation is given by Rosenbluth et aJ. 13 which contains a term of the form 

V-[/ Q (v)VV^(v)]. 

This form is used in several numerical codes, 2-4 even though more derivatives of ipb must be taken. One 
form may be derived from the other by noting that W 2 tp b = 4>b-) 

There is an efficient method for calculating the Rosenbluth potentials. This involves decomposing f b 
in Legendre harmonics eq. (15). Then we have 13 



1 



21 + 1 
1 



2(4Z 2 - 1) 



" fi l \v')dv' 



(25a) 



n'l-3 



l + lv' 2 



(25b) 



Let us assume that f b may be represented by K Legendre harmonics (i.e., the upper limit in the sum in 
eq. (14) is K — 1). Then the calculation of (v) from /b(v) using eq. (15) takes 0(N 2 ) computations 
for each I or 0{KN 2 ) computations altogether. Given f^\v), the calculation of (j) b \v) and ^\v) 
using eqs. (25) takes 0(N) computations for each I. The calculation of </>(v) and V>(v) takes a further 
0(KN 2 ) step. Overall the number of steps is therefore 0{KN 2 ). Often, we can take K to be quite 
small (usually K < 10), and in any case we have K < N, so we can compute the collision term much 
more economically than using the Landau operator directly. 
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IV. Approximations to the Collision Operator 
A. Isotropic background 

If the background distribution is isotropic ,/h(v) = fb{v), then so too are <j> and ip. The collision term 
is then from eqs. (24), (25), and (11) 



ca/6 _ _r>a/b^fa_ _j_ pa/bf 

cv cvv o 1 cv Ja>) 



where 



rya/b _ 

cvv 


47iT a / 6 / 


s: 


3n b \ 




n a/b 
U cBB — 


47iT a / b / 
3n b \ 


j; 


pa/b _ 

cv 


4 7r r°/ b 


m a 


■in h 


mb 



- cvv 1 * cv J a ! (26a) 

S c0 =- D cee--QQi ( 26b ) 



—f b (v')dv'+ v'f b (v')dv'), (27a) 

— (3v 2 -v' 2 )f b (v')dv' + J v'f b {v') dv'j , (27b) 
f v "\v 12 

/ ^n-hW) dv'. (27c) 



na/ 6 


_ -pa/6 ^tfc 

Q ) 

tr 




n a/6 


_ T a/b 1 ( x 

2v \ 


v 2 J 


pa/b 

cv 


= _ Y alb^}_ 

rrib v 2 


1 



B. The high- velocity limit 

If v is much greater than the thermal velocity of particles of species b, the indefinite limits in eqs. (27) 
may be replaced by infinity to give 

-,,2 

(28a) 
(28b) 
(28c) 

where the thermal velocity is defined for an arbitrary isotropic distribution as 

A POO 

< = ^~ v'fs^dv. (29) 
[For a Maxwellian distribution, eq. (18), this reduces to the usual expression eq. (16).] 

C. Maxwellian background 

If the background distribution is Maxwellian eq. (18), the integrals in eq. (27) can be carried out to 
give 15 

/b = ^o_^/b v = _ uerf/(u) n (3()c) 

m a + m b v 2 m b L WJ 
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where 

2 f u 

erf(u) = —= / cxp(— x )dx, 
V 7 *" Jo 

erf'(u) = -^cxp(— m 2 ), 



V2vtb 

The parallel diffusion rate v^ 1 ', perpendicular diffusion rate v± b , and the slowing down diffusion rate 

v^ h are the same as those defined in the NRL Plasma Formulary. 12 [However, the NRL Plasma Formu- 
lary (1983 edition) has an incorrect formula for the collision operator with a Maxwellian background.] 
For u > 0, erf (u) is given approximately by 10 



fc=i 

where 



erf(u) = 1 — cxp(— u 2 ) ^ aut k , 
k 

t= 1/(1 +pu) 



and 

p = 0.32759 11, 

ai = 0.25482 9592, a 2 = -0.28449 6736, 
a 3 = 1.42141 3741, a 4 = -1.45315 2027, 
a 5 = 1.06140 5429. 

This approximation cannot be used in evaluating eqs. (30) near u = because there is cancellation to 
leading order in all three terms. In that case, the Taylor expansion, 



erf(u) - «erf(«) = 4= Y ( ^\, 2 u 2k+3 
w w VtF^ fc 2fc + 3 

v i — n 



^(5» s -!» 5+ r 7 -if" s 

may be used. Alternatively, we can compute eqs. (30) by numerically evaluating the integrals in eqs. (27) 
with fb(v) = fbm{v)- This method is then easily extended to include relativistic effects as given in 
sec. X. 

The ratio FcJ b /Dclv from eqs. (30) [and also (28)] is —m a v/Tb so that the effect of collisions with 
species b is to make species a approach a Maxwellian with temperature T b . 



D. Linearized collision operator 

In many applications in plasma physics (including those involving rf waves) collisions dominate the 
thermal particles. Therefore, the distribution function may be expanded about a Maxwellian 

/o(v) = fam(v) +/ol(v). 

The self-collision operator C(f a , fa) may be approximated by the linearized operator 

Cnn a (/a(v)) = C(f a (y),f am (v))+C(f am (v)J a (v)), (31) 
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where we have made use of the fact that C(f am , fam) = 0, and we have ignored terms of order f^. We 
can compute C(/ a (v), f am (v)) using eqs. (26) and (30). To compute C(/ am (u), / a (v)), we express 
/ a (v) as a sum of Legendre harmonics, eq. (14), to give 



C(f am (v)Ja(v)) = E C (fam(v), /J° W (COS 9)) . 



1=0 



The zeroth term in the sum can be computed using eqs. (26) and (27) giving 

C(f am (v),.i 0) (v)) _ 4rf-/ a 



fam{v) 



+ 



"to 



/ 








v' J 



The next term in the sum in eq. (32) is given by eqs. (24), (25), and (1 1) 



.f am (v)cos9 ~ n a [ L {V)+ J vlJ a {V) 

r°° i/ 2 

Jv v ta 



,,'3 



5v 2 a w 2 3i> 2 



The corresponding fluxes for this term are 

47rr a / a 



r<a/a 
fam (v) COS 



o a / a 

/ am (v)sinfl 



/iV) 



t/ 5 2v' 3 \ 



5w 2 a w 3 3v 3 / 



5v? a v> 2 3w' 2 



dv' 
dv' 



3 6v 2 „v 3u 3 



(32) 



(33) 



(34) 



(35a) 



(35b) 



tion 



are 



Because of conservation of number, momentum, and energy, the solutions to the homogeneous equa- 

^(/ (V))=0 

/o(v) = (a + b • m a v + d\m a v 2 )f am {v). 

If we substitute a = b = 0, we obtain a check on eq. (33). Similarly, a = d = and b = vm gives 
a check on eq. (34). Such checks are useful when incorporating the linearized collision operator into a 
numerical code. 



E. Electron-ion collision operator 

We now turn to the specific problem of current drive by lower hybrid waves. In this problem we wish 
to solve the Fokker-Planck equation for the electrons including the effects of electron-ion and electron- 
electron collisions. 
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Because the ions are so massive relative to the electrons, we have v » vu for nearly all the electrons 
and eqs. (28) apply. Indeed, we can make the further approximations nrn — > oo, vu — > 0, in which case 
the collision operator is given by eq. (26) with 

D e c ii = = 0, (36a) 
^=r e / e g, (36b) 

where 

% lnA e / 4 



q e \nA e / e ' 

and we have assumed neutrality q e n e + q^ = 0. The full electron-ion collision term C(/ e , /j) can be 
written as 

For a multispecies plasma is replaced by Z e g where 

£ s n s g s 2 lnA< 8 

^efi — 



e q 2 e \nk e / e 



and the sum extends over all the ionic species. 

With this collision operator the ions are characterized by a single dimensionless parameter Z t (or 
Z e s). The collision operator allows momentum to be transferred from the electrons to the ions, but there 
is no energy exchange. The non-negative nature of f e is preserved. 



F. Electron-electron collision operator 

There are several choices for the electron-electron collision operator. We will discuss them starting 
at the most complex. 

The full electron-electron collision operator is given by eqs. (24) and (25). This was first used in 
current-drive studies by Harvey et aJ.. 16 This collision operator conserves both momentum and energy. 
The electron distribution f e remains non-negative. Because the collision operator conserves energy, there 
is nowhere for rf energy absorbed by the electrons to go. The problem arises because we have reduced 
a problem in configuration and velocity space to one in velocity space alone, so that there is now no 
spatial diffusion of energy. In practice, this problem is solved by inserting an energy loss term into the 
Fokker-Planck eq. (1). Unfortunately, there are several different models for this loss term and so this 
procedure is somewhat ad hoc. 

The linearized electron-electron collision operator is given by eq. (31). This too conserves momentum 
and energy. The non-negative nature of f e is no longer preserved. When the perturbation to f em is small, 
f e usually only becomes negative far out on the tail. The energy conservation of the collision operator 
again necessitates the introduction of an energy loss term. Fortunately, there is a systematic way to do 
this within the context of a Chapman-Enskog-Braginskii expansion. 17,18 The energy loss term has the 
form 




m e v 2 3\ dhiTg 
-fir-- 2j Um[v) ~dT 



which appears on the right-hand side of eq. (1). Operationally, d\nT e /dt would be adjusted to ensure 
that the energy of the electron distribution / e (v) remained constant. (In fact, one of the results of the 
expansion procedure is an equation for the evolution of T e including the effects of rf and ohmic heating 
and of energy transport.) 
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A useful modification of this collision operator is the truncated collision operator 

C t e / u e nc (/ e (v)) = C(/ e (v), f em (v)) + C(f em (v), fW(v) cos6), (38) 

where the first term is given by eqs. (26) and (30) and the second term by eq. (34). This differs from the 
linearized operator in that we only retain the I = and I = 1 terms in the sum in eq. (32) and we further 
approximate /j ^ (v ) by f em (v). As a consequence, this operator conserves momentum but not energy; so 
there is no need to introduce an energy loss term. Again, the electron distribution function may become 
negative. This collision operator is useful in the study of current drive by low-phase-velocity waves and 
in the treatment of problems with an electric field. In both of these examples, a momentum-conserving 
electron-electron collision operator is required. This operator was used (in a relativistic form) in the study 
of current drive by fast waves. 19 

A slightly different technique for ensuring momentum conservation was used in our study of current 
drive by low-phase-velocity waves. 20 There we approximated the electron-electron collision operator by 



Cdrift(/e(v)) = C(/ e (v),/ ero (|v-^V|||)), 



where the background is a drifting Maxwellian with a drift speed Vd adjusted so that the parallel force 
between f e (v) and the drifting Maxwellian, 



vanishes. This collision operator conserves momentum (by construction) and preserves the non-negative 
nature of f e . Energy is not conserved. It is, however, slightly less accurate than the truncated operator. 
In particular, while the truncated operator gives the correct value for the electrical conductivity, 21 this 
operator gives an answer which is in error by about 15%. The computation of this collision operator 
involves computing eqs. (30) in the drifting frame, converting to cylindrical coordinates using eqs. (10), 
transforming to the rest frame (which is easy in cylindrical coordinates), and finally converting back to 
spherical coordinates using eqs. (10). In order to determine the drift speed, we use the analytical formula 
for the force on an electron Maxwellian drifting with speed v d due to a stationary ion background, i.e., 



u e/i „ 1 /2 / rrn 

P» = -n e m e i/ te v d Zi- 



3 V 7T y nii + m e 

which is valid for \vj\ <C v te - Here we have taken the mass ratio m e /nii to be finite and have assumed 
that Ti = T e . The force between two electron Maxwellians with a relative drift of Vd is found by taking 
Zi = 1 and rrii = m e which gives 

e/e _ n e m e v te v d 
F W ~ 3^ ' ( } 

In the numerical code P^ e is computed using eq. (39). Equation (40) is used to estimate the change in 

Vd required to give P^ e = 0. 

The situation may be further simplified by assuming that the background electrons are Maxwellian, 
so that the collision operator is given by 

Ct(/eW)=C(/e(v),/™W), (41) 

which may be evaluated using eqs. (26) and (30). This operator conserves neither energy nor momentum. 
It does preserve the non-negative nature of f e . It was used by Kulsrud et al. 5 in the study of runaways, 
and in studies of lower hybrid current drive. 6 The Maxwellian background serves as a heat bath, so 
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no energy loss terms are required. The absence of momentum conservation introduces approximately a 
factor-of-two error in the electrical conductivity 5 and in the efficiency of current drive by slow waves. 20 
There is a relative error of order (v te /v) 3 in the determination of the current-drive efficiency for fast 
waves. 19 

Lastly, may be approximated by using the high velocity limit, i.e., by using eqs. (28) instead of 
eqs. (30). In fact, because eq. (28b) gives negative diffusion for small v, it is usually replaced by 

n a/b _ -pa/6 1 

We define the resulting electron-electron collision operator as C^ h . It has much the same properties as 

^Max- 1° particular, it yields a Maxwellian (with temperature T e ) as the steady-state solution. Because 
of the greater error in the collision term for thermal particles the electrical conductivity is even lower 
than for C^ x . The evaluation of eqs. (28) is, of course, a little easier to program than that of eqs. (30). 
However, because the results of evaluating eqs. (30) can be stored in a table, the extra computational cost 
of working with C^tx ls insignificant compared to the solution of the Fokker-Planck equation. Since 
C'high is less accurate, its use is not recommended for numerical work. It is, however, useful in analytical 
work. 

When working with these electron-electron collision operators, it is useful to have some benchmark 
against which to check their numerical realization. A useful benchmark is provided by the electrical 
conductivity, which is the ratio of electrical current to electric field in the limit E — > 0. This is tabulated in 
table I for various values of Zi and for all the electron-electron collision operators discussed here. These 
values were obtained by solving the corresponding one-dimensional equation by the method outlined in 
sec. XI. In all cases, the electron-ion collision operator is given by eq. (37). The conductivity using the 
full and the truncated electron-electron collision operators is the same as for the linearized operator. In 
the limit Zi — ► oo, the conductivity is independent of the electron-electron collision model 

J [2 1 n e ql 



E V 7T Zi m e v te ' 

For the high-velocity approximation to the collision operator the conductivity can be expressed analyti- 
cally as 

J 1Q [2 3Z, t + 13 n e ql 



E 3 V n (Zi + 3)(Zi + 5) m e v te 



V. Quasilinear Operator 
A. Single wave 

The interaction of electrons (or other species) with a wave is conveniently described in terms of the 
quasilinear theory. 22 In this theory the flux of electrons in velocity space is given by 

S w = -D w -Wf e , (42) 

where is the quasilinear diffusion tensor which depends on the waves present in the plasma. Although 
quasilinear theory is not strictly applicable to a single wave, we will start with this case because it is the 
simplest. Suppose there is a uniform wave present in the plasma, i.e., 

E(r, t) = RcfE^ exp(ik • r - iut)} . (43) 
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The quasilinear diffusion coefficient is given by 

7T q 2 
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and 



Dw = o ~~% <y ( w ~ k W v W ~ n ^K a ™ (44) 



fc|| 

U) 



V\\ ]V_L + Dj_V|| 



ft _ E w+ J n -i + E w -J n+ \ ^|| 

V2 d_l 11 

where £l e = q e B/m e is the electron cyclotron frequency, £> is the magnetic field, * indicates complex 
conjugation, J n is the nth order Bessel function, and the argument of the Bessel functions is kj_v±/Fl e . 
E w+ and E w - are the left- and right-handed components of E^; in a right-handed cartesian coordinate 
system with z parallel to B and k lying in the (x, z) plane, we have 



E w + 



E„ 



Ewx H~ 



ivy 



V2 



E — i E 



V2 ' 
E w \\ — E wz . 

It is instructive to consider the properties of eq. (44). The delta function specifies the resonance 
condition. Only particles for which the Doppler-shifted wave frequency w — kuVu is zero (n = — the 
Landau resonance) or a multiple of the cyclotron frequency (n ^ — a cyclotron harmonic resonance) 
interact with the wave. The vector a„ is perpendicular to the velocity of the electron in the wave frame 
v — (w/fc||)v||. This means that the wave-induced flux is along diffusion paths which lie in constant- 
energy surfaces in the wave frame; see fig. 2. Similarly, the flux is proportional to the gradient in f e in 
this direction. As a consequence, when an electron interacts with a particle via the Landau resonance, the 
diffusion tensor consists of only a single component 

D t0 =D||||V||V||. 

Likewise, for a cyclotron harmonic resonance, we have 

D w = Duvivi, 

provided that v± is small compared with nQ, e /k\\ . 

The appearance of the delta function in eq. (44) is a consequence of the assumed uniformity of the 
magnetic field. In this case, v\\ is a constant of the unperturbed motion and so a particle remains in reso- 
nance for a long time. In situations described by bounce-averaged codes, the magnetic field and dm vary 
along a particle orbit so that the particle does not remain in resonance. This effect can be taken into ac- 
count by averaging eq. (44) along a particle trajectory. 23 This removes the delta function, although there 
are still singularities in the resulting expression arising from those particles which turn in the resonance. 8 



B. Many waves 

Equation (44) is easily generalized to include a more realistic representation of the wave fields. An 
important application is to the incorporation of quasilinear effects into a ray-tracing code. Here the 
externally injected rf power is represented by several rays. Let us consider the interaction of these waves 
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with the electrons on a particular flux surface. At the point where a given ray intersects the flux surface 
it is characterized by its position r, wave number k, and power W (usually the frequency w is fixed by 
the rf source). W measures the number of watts carried by the ray. In order to apply eq. (44), we must 
determine the amplitude of the corresponding single wave E w which has the same polarization and same 
rms field amplitude as the ray (with the rms averaging performed over time and over the flux surface). 
The ray contributes 

W 

U=- — - 

W 9 ■ n| A f 

to the wave energy density (in J /m 3 ) averaged over the flux surface, where v g is the group velocity of 
the ray, n is the unit vector normal to surface at the point of intersection and A / is the area of the flux 
surface. The polarization of the electric field is given by 

K • E w = 0, 

where 

K = 4(kk-fc 2 l) + l + ^^ 

is the dispersion tensor, c is the velocity of light, and cr(k, uo) is the conductivity tensor. The energy 
density is related to ~E W by 24 

Given W, we can therefore determine E w (to within an ignorable phase factor) appropriately averaged 
over the flux surface. 

This is substituted into eq. (44) and the result summed over all the rays to give the overall quasi- 
linear diffusion tensor. In practice, the delta-functions appearing in this expression must be replaced 
by smoothed functions. This allows the ray-tracing procedure to reflect the true situation in which a 
continuous spectrum of waves is launched. 

We complete the discussion of the ray-tracing by pointing out that the damping of the rays should 
be calculated self-consistently from D w . The power that a particular ray loses per unit volume due to 
absorption by the electrons is given by eq. (20), where instead of the total we use the contribution 
the ray in question makes to S w . To this should be added the power absorbed by the other species if 
applicable. Then the ray power W satisfies the equation 

dW 

- w = -P\v B -n\A„ 

where the time derivative is the derivative taken along the ray. 

If instead of a discrete set of waves, the wave fields are given by a spectrum 

, ... . . , d 3 k 

E 



3(r, t) = J E w (k) exp[ik • r - iuj(k)t] 



(2tt) 3 ' 
then eq. (44) becomes 22 

v— n Q 2 f rf 3 k 1 

D - = E 4 J p^^M k ) - k U v \\ - nfleKa™, (45) 

where V p is the configuration space volume of the plasma and the definition of a„ is generalized in the 
obvious way. 
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C. Model forms 



The results given above allow a ray-tracing code to be coupled to the solution of the Fokker-Planck 
equation. This is an extremely complicated system, and much work has been carried out using assumed 
forms for the quasilinear diffusion coefficient. This allows us to study the physics of the interaction of the 
electrons and the waves without having to worry about the additional physics of the wave propagation. 
The most widely used model form for lower hybrid waves was introduced by Fisch 1 and is given by 

O w = A„(v||)v||V||, (46a) 

where 

This form of D„, is justified as follows. Because lower hybrid waves interact only via the Landau 
resonance, only the V||V|| component is present. If k_iv te /il e <C 1, the dependence on perpendicular 
velocity may be ignored (J ~ 1). Finally, in many cases, the magnitude of the quasilinear diffusion 
greatly dominates over the collisions; thus the quasilinear diffusion coefficient tends to make an abrupt 
transition (in velocity space) from being negligible to being large; if D is sufficiently large (i.e., large 
enough to form a quasilinear plateau), this situation is accurately modeled by eq. (46b). 

This particular form for D^, is useful because much theoretical work has been carried out using it. 1 
Numerical solutions to the Fokker-Planck equation provide the best test of these theories. It is therefore 
important that any numerical code be able to handle the discontinuities in D w . (Note, however, that both 
/ e and S are continuous even if is not.) 

This model is readily generalized, for example, by allowing Dq(v\\ ) to be an arbitrary function. Thus 
the effect of a backward component to the lower hybrid spectrum can be studied by including another 
boxlike component to D . Similar models have been used to study low-phase-velocity current drive 20 
and electron-cyclotron current drive. 25 



D. Direct specification of the quasilinear flux 

Both analytical and numerical studies show that the current drive efficiency is primarily determined 
by the location at which electrons interact with the waves and the direction in which the waves push the 
electrons. It is sometimes useful to specify the rf-induced flux directly as some arbitrary vector field 
S w (v). Indeed, in some cases we may know more accurately than we know D^,. In a ray-tracing 
calculation, may be calculated self-consistently in terms of the power flows in the various rays. 
However, in cruder zero-dimensional calculations, we may wish to assert merely that so much rf power 
is absorbed by the electrons. Then S„, may be estimated from eq. (20) using an a priori knowledge of 
which electrons interact with the waves. Alternatively, S w may be estimated from either an approximate 
analytic solution of the Fokker-Planck equation 26 or from a solution of the one-dimensional Fokker- 
Planck equation. 1 

If S w is given, then the Fokker-Planck equation (1) is an inhomogeneous (instead of homogeneous) 
equation. However, assuming that one of the linear electron-electron collision operators is being used, 
the linear operator acting on f e in eq. (1) is now independent of the wave drive. This property is used in 
the adjoint methods to provide a very efficient method of solving for moments of f e (see sec. XI). 
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VI. Boundary Conditions 



A. Computational domain 

We shall take the computational domain V for the Fokker-Planck equation to be 

< v± < v± max , U|| min < W|| < W|| max , (47) 
for problems solved in a cylindrical coordinate system and 

< v < v max , < 9 < 7T, (48) 

for problems solved in a spherical coordinate system. The boundary of V is defined to be A. (For 
example, in spherical coordinates, A is the spherical surface v = v max .) 



B. Internal boundaries 

We distinguish two types of boundary: internal and external boundaries. The internal boundaries are 
the simplest. In a cylindrical coordinate system we have an internal boundary at v± — 0. Values of / 
beyond this boundary are determined by symmetry 

fe{-V±,V\\) = f e {v±,V\l). (49) 

Similarly, in spherical coordinates we have internal boundaries at v = and at 9 = and 9 = tt. These 
boundaries are treated with the boundary conditions 

f e (-v,6) = f e (v,w-6), (50a) 
f e (v,-6) = f e {v,0), (50b) 
f e (v,7T + 9) = f e (v,7r-9). (50c) 



C. External boundaries 

The other boundaries are inserted into the problem in violation of the true physical picture. In reality 
the velocity domain extends off to infinity; on the computer, however, we normally study only a subspace. 
We have to choose the subspace to include all the interesting physics: for studies of electron distribution 
in a spherical coordinate system, we require v max » v te ; if the electrons are driven by lower hybrid 
waves, then we further require t> max > (w/fc||) ma x, the maximum wave phase velocity; if we wish 
to study runaways, then v max must exceed the runaway velocity; and so on. We next have to choose 
boundary conditions which are as "innocuous" as possible; i.e., which perturb the solution in the domain 
of integration as little as possible compared to the solution in the full domain. 

For electron current-drive problems we choose the condition 

S-n = 0, (51) 

on the external boundary A, where ft is the normal to A. This means that plasma cannot enter or leave 
the domain of integration. Thus the number of electrons is conserved with this boundary condition. This 
boundary condition gives a Maxwellian steady state in the absence of the rf, and allows a steady-state 
solution to be reached in the presence of rf. 

If an electric field is present, then in the real problem some electrons will run away. Now we wish 
to impose boundary conditions which "allow" this to happen. At the boundary we have v » v te so that 
collisions are weak, and the dominant process is the acceleration by the electric field (we assume that the 
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boundary is removed from the region where the rf diffusion takes place). The Fokker-Planck equation 
then reduces to a hyperbolic equation. The tactic is to apply the same boundary condition as before, 
namely eq. (51), where the characteristics of the hyperbolic system enter the domain of integration. 
Where the characteristics leave, we set those diffusion terms which lead to a flux across the boundary 
to zero. This makes the equation purely hyperbolic in the direction normal to the boundary and so no 
boundary condition is required. (We shall see in sec. VII how this comes about in the numerical scheme.) 

If we assume that q e E > so that electrons run away in the positive direction, then in cylindrical 
coordinates we would impose 

S|l = 0, fort>|| = U|| min , 

S± = 0, forv± = v ±max , (52) 

D ll-L = ^11 II = °. forw [| = U ||max- 

(The boundary at v± = v± niax is taken to be an incoming boundary because the small collisional friction 
makes the characteristics enter along this boundary.) 

A slightly more accurate treatment is possible in spherical coordinates. If we compare the various 
collision terms in the high-velocity limit eqs. (28), we find F cv <~ D c gg/v ~ 1/v 2 and D cvv /v <~ vf e /v 4 . 
Thus we can ignore the energy diffusion term D cvv compared with the other collisional terms. The 
pitch-angle scattering term D c gg requires no special handling because it causes diffusion parallel to the 
boundary. The equation is, therefore, hyperbolic in the direction perpendicular to the boundary with a 
characteristic acceleration given by F v = F cv + (q e E/m e ) cos 9. The boundary conditions on v = w max 
then become 

S v = 0, for F v < 0, 

(53) 

D vv = D v g = 0, for F v > 0. 

For w max » v te , F cv is accurately approximated by eq. (28c) (with a = b = e). Thus, if \E\ < 
m e r e / e / \q e \ u^ax' this boundary condition reduces to eq. (51), allowing problems involving both an 
electric field and rf diffusion to be handled in a unified way. In this small electric field limit, S v = is 
zero everywhere on the boundary and the numerical runaway rate vanishes. This is a close approximation 
to the true situation in which the runaway rate is exponentially small — on the order of exp(— v^^/v^). 



D. Treatment of runaways 

With a finite boundary, we can determine the runaway rate accurately (provided w max is sufficiently 
large). However, the behavior of the runaways beyond the boundary is not followed. One could, of 
course, just choose a very large boundary; but this is wasteful of computer resources and really just 
postpones the time at which the problem is encountered. It is, therefore, preferable to treat the runaways 
as a separate species. Assuming that the runaways are affected only by the electric field, the density and 
current moments of the runaway population form a closed set of equations. We define 




where V is the complement of V, i.e., the region v > v max in spherical coordinates. Applying eqs. (5a) 
and (5b) to V we find 

= ^n r + / q e v l{ S-d 2 A. 
dt m e J A 
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Thus if we wish to determine the total current as a function of time, we need only supplement the Fokker- 
Planck equation by two ordinary differential equations and then sum the nonrunaway and runaway con- 
tributions to the current. 

VII. Spatial Differencing 

A. Choice of coordinate system 

We have discussed both the cylindrical and the spherical coordinate systems. Which one should be 
used in a given application? The numerical scheme that is described here works best if the diffusion 
tensor is nearly diagonal. Then the mixed derivative terms in eqs. (8) or (9) are small. (It is these 
terms which tend to make the numerical scheme unstable.) Now the collision operator is approximately 
diagonal in spherical coordinates while the quasilinear term is nearly diagonal in cylindrical coordinates. 
Thus the choice of coordinate system to some extent depends on the relative strength of these two terms. 
Cylindrical coordinates were used in the study of current drive by low-phase-velocity waves 20 because 
the edges of the resonant region line up with coordinate lines allowing the scaling with phase velocity 
to be measured more accurately. On the whole, however, the spherical system is to be preferred because 
the electron-ion collision term eq. (37) becomes large near v — and we wish this term to be diagonal. 
In ref. 20 much smaller time steps had to be taken to avoid the problem with the electron-ion term. The 
boundary conditions can also be applied more accurately in spherical coordinates when an electric field 
is present [eqs. (53)]. For this reason, we will focus on the spherical coordinate system in this section. 
Extension to the cylindrical coordinate system is straightforward. 

An alternate representation of f e is as a series of Legendre harmonics. This has no particular merit 
in quasilinear problems because the sharp gradients in D„, eq. (46), cause the Legendre expansion to be 
slowly convergent. 

B. Normalizations 

In solving equations of physical significance on the computer, it is often useful to normalize all 
the physical quantities. This allows us to work with numbers which are closer to unity (and thus avoid 
potential problems due to arithmetic overflow or underflow); more importantly, the number of parameters 
needed to specify the problem is often reduced. 

For the problem of current drive by lower hybrid waves, we solve the Fokker-Planck equation for 
the electrons. We normalize velocities to Vt e eq. (16), times to T te eq. (17), the electron density to n e , 
the electron distribution to n e /vf e , the quasilinear diffusion coefficient to v1 e v te , the electric field to 
m e vteVtel 'q e , the current density to n e q e vt e , power density to n e m e v1 e vt e , etc. 

These normalizations coincide with those used by Kulsrud et ah. 5 However, they differ from those 
used in some of our earlier papers, e.g., ref. 6. (The thermal collision time differs by a factor of two.) 

Since we are only dealing with the electron distribution, we will drop the species label from / and 
other electron quantities. Otherwise, we shall use the same notation for normalized and unnormalized 
quantities. For example, the electron Maxwellian eq. (18) reads in normalized terms 

= ^7jexp(-ii; 2 ). 

The reduction in the number of parameters now becomes apparent. The plasma is characterized by a 
single parameter Zj and the quasilinear diffusion coefficient by three parameters D , V\, and v 2 - 
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C. The numerical grid 

We wish to solve eq. (1) in the domain V eq. (48). We do this by converting the differential equation 
to an algebraic equation using the finite difference method. In this method / is represented by its values 
on finite set of points and differentials are represented by differences between neighboring values. 

First, we establish a numerical grid by dividing v and 9 into N and M equal pieces, respectively. 
Thus we define 

Av = iw/iV, A6 = tt/M, (54) 

together with grid positions 

Vj = j Av, (55a) 
9, = i A9. (55b) 

This grid system defines a system of cells. The electron distribution function is represented by its values 
at the centers of these cells, i.e., by the values 

/i+i/2j+i/2 = /0>j+i/2. 0i+i/2), ^ < i < M, < j < iV, 

with i and j being integers; see fig. 3. The cell Vj < v < Vj + \, 9i < 9 < 9 i+ i (i and j integers) is 
assigned a volume 

^+1/2^+1/2 = 27rsin0 i+1/2 u? +1/2 AvAO. (56) 
We will define numerical volume integration by 

M-l N-1 

intX= ^2 X] X i+l/2.J + l/2.fi+l/2.J + l/2V i+ l/2,j + l/2- (57) 
i=0 j=0 

This is the discrete analogue of J v Xfd 3 v; see eq. (13). We define the flux of a quantity through the 
boundary by 

M-l 

fluxX = X 2TTsm9 i+1/2 v 2 N X l+1/2 M+i/2S V a + i/2,N AO, (58) 

i=0 

which is a discrete analogue of J A XS ■ d 2 A. The number density of electrons becomes 

n = intl. (59) 



An alternative approach to finite differences is provided by the finite-element method where the / is 
represented by the superposition of a set of trial functions with finite support. This approach has been used 
in Fokker-Planck codes by workers at Lausanne. 27 ' 28 The finite-element method is also used in some 
commercial computer codes for the solution of partial differential equations. One such code has been 
applied to the Fokker-Planck equation by Fuchs et aJ.. 29 If we identify the weights of the trial functions 
with the values of / at the grid positions, we see that the finite-difference and finite-element methods 
are quite similar. In particular, the goals of the methods are identical: to express algebraically df/dt 
at a particular location in terms of / at the same and neighboring locations (usually, the eight nearest 
neighbors). Thus our discussion of the time advancement of the equation in sec. VIII is independent of 
the choice of method. 
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D. Divergence of flux 



Consider the Fokker-Planck equation in the form eq. (2). This is translated onto our numerical grid 
in a conservative form as 

dfi+l/2,j+l/2 _ ( Vj +1 S v s +1 / 2 ,;j + l - V?S v>i+1 / 2 ,j 



dt V " v j+i/2 ^ v 

_l_ sin^ i+ i»S , e ^ + ij + i/2 - sm6 i Sg tit j + i/2\ 
vj+i/2 sin 6 i+1/2 A6 J' 

Notice that the fluxes are required on the edges of the cells (see fig. 3) and that the fluxes on the internal 
boundaries do not contribute since they are multiplied by u =0 or sin O = sin 9 m = 0. With this 
method we difference the fluxes and not the diffusion and friction coefficients. This lets us treat problems 
in which D„, is discontinuous, e.g., as given by eqs. (46). The scheme in eq. (60) is accurate to second 
order in Av and A9. 

This form of difference equation is called conservative because it obeys the conservation law 

^ + flux 1 = 0> (61) 
at 

where int and flux are defined by eqs. (57) and (58). This is a discrete counterpart of eq. (5a). If 
S v ,i+i/2,N — for all i, then we have flux 1 = and particles are exactly conserved in the numerical 
scheme (if we ignore round-off errors). The discrete form of the parallel component of the momentum 



conservation law eq. (5b) is 

— h flux(« cos 9) = ^2 E 27rw i sin ^+1/2 cos ^+1/25^+1/2 ,j AvA9 



M-l N 



i=0 j=0 
M N-l 

_ E E 27 ™f+i/2 sm%S etiJ+1/2 Av 2 sin(i A6), (62) 

2 = j=0 

while the energy conservation relation eq. (5c) becomes 

dmt(±v 2 ) M_1 N 

+ flux(i« 2 ) = E 2n sin ^+1/2^^^+1/2,, AvA9. (63) 

These relations are useful in that they establish definitions of various physical quantities that are consis- 
tent with the numerical scheme. For example, we can interpret the right-hand side of eq. (63) as the total 
power flowing into the electrons. This definition is consistent with the numerical definition of the energy 
of the electrons, namely int(^v 2 ). Furthermore, we can determine the power flowing into the electrons 
from the waves (for example) by replacing S v in the right-hand side of this equation by the flux due to the 
waves S wv [compare with eq. (20)]. In this way, we obtain a complete and accurate power balance for 
the electrons. Similarly, the right-hand side of eq. (62) gives the definition of the force on the electrons. 
This is used when evaluating P'^ e in eq. (39). 

[In order to prove eqs. (62) and (63), the following relation is useful: 

M-l M-l 



E 



+ Ai)(B i+ i — Bi) — AmBm — AqBq — \ 2 (Bi+i + Bi)(A i+ i — Ai). 



i=0 4=0 



This is the rule for "summing by parts" — the discrete counterpart of integration by parts.] 

The basic difference equation (60) is readily generalized to nonuniform grids. However, the deriva- 
tion of eqs. (62) and (63) relies on the uniformity of the grid and they cannot easily be generalized. 
Nonuniform spacing is used in FPPAC. 4 
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E. Stream function 



A very useful tool for understanding the Fokker-Planck equation (2) is the flux plot, which shows the 
vector field S(v). This is sometimes displayed as a set of arrows, one at each grid point, which point in 
the direction of S and which have a length proportional to S. In this problem, S v and Sg are known at 
different locations, so that realization of this prescription would necessitate interpolation. Furthermore, 
such a display is often very misleading because the visual impression is strongly affected by whether the 
arrows line up with other grid points or not — a purely artificial aspect of the problem. 

The much superior method is possible if we restrict ourselves to the steady state. In this case, the 
vector field S(v) is divergence-free V • S = 0, and so may be expressed as the curl of a stream function, 
i.e., 

S(v)=Vx^4, 
2ttv sin 

where <p is the azimuthal coordinate. The components of S are given by 

n dA 
2nv 2 smo o0 
n dA 

Se = — —, ir . (64b) 

Zirv sin v ov 

Because S • VA = 0, lines of constant A are stream lines. Thus a contour plot of A(v) gives the vector 
field of S(v). The stream lines are obviously closed (indicating that the flow is divergence-free), and the 
total flux of electrons between any two contours is equal to the difference in the values of nA on those 
two contours. 

We can compute A on the numerical grid using discrete analogs of eqs. (64) 

2™? 



Aij = 3 - 8in0i, +1 / 2 S Vti , +1/ 2,j AO, (65a) 

n i'=0 

2tt sin 6. 



n 



^ v jl+1 /2Se iifjl+1 /2 Av. (65b) 

j'=0 



If dfi+i/2,j+i/2/dt = according to eq. (60), then these two definitions are consistent. 



F. Computation of the flux 

In order to complete the specification of the difference scheme we must give formulas for S vA+1 / 2 ^ 
and Ses j+1/2 m ec l- (60). These depend on the type of electron-electron collision operator used. We start 
with collisions off a Maxwellian background C^ x , eq. (41). This is the simplest case and yet it exhibits 
all the difficulties of solving the Fokker-Planck equation. 

The collisional flux is given by the sum of the flux contributing to C^ax which is given by eqs. (26) 
and (30) and the flux contributing to C e l % which is given by eqs. (26) and (36). [In fact, we compute 
the electron-electron flux by numerically evaluating the integrals in eqs. (27).] To this is added the 
quasilinear flux from eqs. (42) and (46) and the electric-field-induced flux from eq. (4). Both these 
terms are converted into spherical coordinates using eqs. (10). The total flux is then given by the general 
equations (9b) and (9c). 

The diffusion and friction coefficients are computed at the points at which we need to know S v and 
S . Thus we compute D vv ^ +1 / 2 ,j, A,0,i+i/2,j» F v ,%+i/2,j, and ^e«,»,j+i/2. D$o,i,j+i/2, Fe,i,j+i/2- 
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The coefficients for S v are not required at j = 0, nor those for Se at i = 0, M, because these fluxes are 
multiplied by zero in eq. (60). The boundary conditions eqs. (53) at w max are handled by setting 

Dyv,i+1/2,N 0, 

D v $,i+\/2,N <— 0, 

Fva+1/2,N <~ max (^,i+l/2,Af > 0). 

Next we must specify the way in which / and its derivatives are to be computed at the edges of the 
cells — i.e., locations (i + 1/2, j) and (i, j + 1/2) — in terms of the values of / at the centers of the cells 
(i + 1/2, j + 1/2). Two of the terms are straightforward: 

dfi+l/2,j /»+l/2,j+l/2 - /i+l/2,j-l/2 

, (66a) 



dv Av 

d/i J + l/2 _ /i+l/2,j+l/2 - /i-l/2,j+l/2 

oe ~ AO 



(66b) 



and 



Again these expressions are accurate to second order. 

The evaluation of / at the cell edges uses a method proposed by Chang and Cooper 30 extended here 
to two dimensions. The simple method, i.e., 

fi+l/2,j = 2(/*+l/2,j+l/2 + /i+l/2,j-l/2)) 

turns out to give poor results for the steady-state distribution. Chang and Cooper replace this with 

/i+l/2J = (1 — ^i+l/2,j)/i+l/2,j+l/2 + $i+l/2,jfi+l/2,j-l/2i (67a) 
/ij + 1/2 = (1 — ^j + l/2)/j + l/2j + l/2 + ^i,j+l/2/i-l/2,j+l/2) (67b) 

where the 5s are given by 

t>i+i/2j = g(~AvF v ^ i+1 /2j/D VVti+1 / 2 .j), (68a) 
<*i,j+i/2 = 9(-A8Fo,i,j+i/2/D e e,i,j+i/2), (68b) 

= tV^- (69) 

w cxp(w) — 1 

The role of the S is to weight the averaging performed in eqs. (67). The weighting is needed because 
often / is a strongly (exponentially) varying function of v. An acute example of this is the Maxwellian 
distribution which varies very strongly for large v. In fact, the weighting is such that a Maxwellian is an 
exact steady-state solution when there is no rf and no electric field and when C^ x is employed as the 
electron-electron collision operator. This is easily seen because for any isotropic distribution S c g = 0; in 
that case, we also require S cv = in the steady state (because there are no sources or sinks of electrons). 
Using eqs. (26a) (with a = b = e), (66a), and (67a), together with F^ e i+1/2 j/D e J v e v i+1/2 . = -vj, we 
find 

/*+l/2,j+l/2 /m,j+l/2 . a x 

7 — ' — = - f = exp(-u j Av). 

Ji+l/2,j-l/2 Jmj-1/2 

The errors in various moments of / are, therefore, exponentially small. With one-dimensional equations 
the weighting cures the problem of / becoming negative. 30 With our two-dimensional equation, this 
problem is alleviated but not cured. In general, this problem is solved by taking a sufficiently fine mesh 
(assuming that the electron-electron collision operator preserves the non-negative nature of /). 
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The function g has the properties 



g(w) = l-g(-w), 

, . 1 w w 3 
• 9H= 2 - 12 + 720 + --" 

ff(-oo) = 1, «?(0) = ^ 0(00) = 0. 

The first two properties are useful for evaluating g(w) for w ^ 1 and ujrjO, respectively. 

The values of the cross-derivative terms which multiply the off-diagonal terms in the diffusion tensor 
(D v g and Dg v ) are now given in terms of eqs. (67) as 

d fi+l/2,j _ /i+3/2j - fi- 1/2 ,j 



89 2A6 

dfij + 1/2 _ fi,j+3/2 - fij-1/2 

dv ~ " 2Av 



(70b) 



The internal boundary conditions eqs. (50) give the values of fi+i/ 2 .j+i/ 2 beyond the internal bound- 
aries as 

/i+l/2,-1/2 = /m-i-1/2,1/2; 
/-1/2J + 1/2 = /l/2,j+l/2) 
/m+1/2j + 1/2 = /m-1/2j' + 1/2- 

These conditions are only needed in the evaluation of cross-derivative terms. The form of eq. (60) auto- 
matically takes care of the internal boundaries for the other terms. 

The external boundary at v = w max is treated as follows: In the computation of SV>,i+i/2.jv we nee d 
only worry about the friction term (since D vv = D v9 = on the boundary) so that only /j+i/2,jv 
is needed. Furthermore, the friction coefficient F v ^ + i/ 2 .n is non-negative. From eq. (67a), we have 
/i+i/2,jv = fi+i/2,N-i/2 because S i+1/2 ,N -» 1 for F vA+1/2 ,n > and D vvA+1/2 .n = 0+. (Obvi- 
ously the value of fi+\/ 2 .N is not required where F v _ i+ i/ 2tN = 0.) Recall that the equation reduces to 
hyperbolic type on this boundary, so that no boundary condition should need to be specified here, as in- 
deed is the case. In fact, the method reduces to the standard upstream differencing for a hyperbolic equa- 
tion on this boundary. In the computation of Sg titN _i/2, only the cross-derivative term df itN _i/ 2 /dv 
potentially involves points outside the integration domain. In this term, we use 

9fi,N-l/2 _ fi,N-l/2 — fi,N-3/2 

dv Av 

instead of eq. (70b). 



G. Matrix formulation 

For collisions off a Maxwellian background the problem is linear so that eq. (60) can be rewritten as 

%+Af = h, (71) 

where / is a vector of length MN of the values j "i+i/2,j+i/2 an d A is an MN x MN matrix of coef- 
ficients. The right-hand side h (also a vector of length MN) is inserted to aid in the treatment of other 
collision operators. For the Maxwellian collision operator, we have h — 0. It is convenient to split A into 
three pieces, namely 

A = A V + A e + A x , 
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where A v contains the terms proportional to D vv and F v , Ae contains those proportional to Dee and F$, 
and A x contains the cross-derivative terms proportional to D v $ and D$ v . With the difference scheme 
given in this section A v and Ae are tridiagonal matrices. Thus we can write 

i+l/2,j+l/2 - a v ,i+l/2,j+l/2fi+l/2,j-l/2 + K,i+l/2,j+l/2fi+l/2,j+l/2 

+ Cw,i+l/2,j+l/2/i+l/2,j+3/2! (72a) 
{■A#f)i+l/2,j+l/2 = tte,*+l/2,j+l/2/i-l/2,j+l/2 + &9,t+l/2,j+l/2/i+l/2,j+l/2 

+ c e,i+l/2,j+l/2/i+3/2,.j + l/2i (72b) 



where 



U J /" ^vv,i+l/2,j 



a-v,i+i/2,j+i/2 — -g- ^ ^ — ! ^,i+i/2,j^i+i/2j ) , (73a) 

_ V 3 ( AnM+l/2J „ 

0w,i+l/2,j+l/2 - "g-l ^ ^Vi+l/2,.j e i+l/2,j 

+ "#7 ( ^'^ /2 ' 3+1 + F v , l+1/ 2, J+ i&. 2,,,-u.i j (73b) 

<Vi+l/2j + l/2 - -g— I ^ H *t),i+l/2,j+iei+l/2,j+l j i (/JC) 

_ sin 6^ / Dee,i,j+i/2 ,. , > 

ae,»+i/2,j+i/2 - -g^- ^ — W . +1/2 A6» ~ ^«.*.i+ 1 /2 () »,j+i/2 ) > ( /3d J 

_ sin #j / 

0e,i+l/2,j+l/2 - -5— — TV - ^9,i,j+l/2ei,j+l/2 



+ (, v j+1/2 A0 + ^^+i^+V2^ + i/2 ) ■ (73e) 



_ sin0 i+ i / D e g ti+1>j+1 / 2 \ 

Ce,i+l/2,j+l/2 - g^ ^ Vj + 1/2 A6 ^0,i+l,j+l/2^i,j+l/2 J , (/3t) 

where e = 1 — 5, B v = Avv? +1 , 2 , and Be = Vj +1 / 2 A9 sin# i+1 / 2 - With these coefficients the bound- 
ary conditions are reflected in the relations 0^+1/2,1/2 = c v,i+i/2,N-i/2 = and 0^1/2^+1/2 = 
c e,Ar-i/2j+i/2 = 0, which are automatically satisfied. 

The matrix A x is more complicated with (^4 x /)j+i/2,j+i/2 depending, in general, on the eight 
nearest neighbors to /i+i/2,j+i/2- The boundary conditions have to be explicitly included in this matrix. 
We do not give expressions for the components of A x here because only the product A x f is ever needed 
in the calculation. This is most easily computed directly in terms of the flux; this also cuts down on the 
storage requirements. 



H. Alternate collision operators 

The methods we will describe in the next sections for solving eq. (71) depend on the linearity of this 
equation and the fact that A v and A e are tridiagonal matrices. With more complicated electron-electron 
collision operators, these conditions no longer hold. However, the techniques can still be used because 
the difference between the other collision terms and the Maxwellian collision term varies slowly in time. 

If the full electron-electron collision operator is used, the basic framework given above still applies, 
except that the diffusion and friction coefficients D e J e and ~F e J e are now given in terms of gradients of 
the Rosenbluth potentials eqs. (24). These coefficients depend on / making the equation nonlinear. In 
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practice, the dependence on / is weak so that the coefficients only need to be recomputed occasionally. 
This also means that the equation is approximately linear so that the linear matrix techniques used to 
advance the equation in time still apply. 

If the linearized or truncated collision operators are used, then the equation remains linear but with 
a term which involves an integral over /, namely C(/ m (w), /(v)) or the truncation of this term. Again, 
this term is weakly dependent on / so that it need not be recomputed every time step. It is then most 
convenient to regard this term as the inhomogeneous driving term h eq. (71). For the truncated col- 
lision operator C t e / U e nc , eq. (38), the elements of h are given by C(f m (v)J^(v) cos 6) evaluated at 
( v j+i/2, 81+1/2)- The computation of this term is described in appendix A. 

VIII. Time Differencing 

A. Crank-Nicholson method 

We now turn to the method for advancing the Fokker-Planck equation in time. If the time step is At, 
then we define 

f k = f(t = t k ), t k = kAt. (74) 
The simplest way of advancing eq. (71) is the explicit scheme 

fk+1 _ fk 

3 - — —L- + Af k = h. 
At 

This is only accurate to first order in At. Furthermore, At must be chosen to be very small, on the order 
of Av 2 or AO 2 , for stability. These defects are easily remedied by the Crank-Nicholson scheme 31 which 
reads 

fk+1 _ fk fk+1 1 fk 

This scheme is accurate to second order in At and is stable if A is positive definite. (This is a condition 
possessed by the continuous form of the operator A.) In order to solve eq. (75) for f k+1 we have to 
compute the inverse of (7 + i At A). This is a large banded matrix which can either be inverted using 
iterative methods or using Gaussian elimination. In both cases the number of operations is 0(N 3 ), 
(assuming M ~ N) making it a very expensive proposition. (This approach is discussed further in 
sec. IX.) 



B. Alternating-direction-implicit method 

Although (I + \At A) is difficult to invert, the matrices (I + \ At A v ) and (I + \At A e ) are rather 
easily inverted. This allows the alternating-direction-implicit method 31 to be used. Unfortunately, (I + 
\ At A x ) is not easily inverted and this means that the cross-derivative terms are treated explicitly in this 
method. Consider the equation 

I + ^A v ) ( I + ^Ae] -^—^ + A.f = h. (76) 



2 7 v 2 / At 

If we rearrange the terms in this equation to give 



^ 2 , „ \f k+1 -.f k , fA , A J k+1 +f k , ., fk 



I + — A v A e J J — - + (A v + A e )- ^-i- + A x f = h, 

we see that this method differs from the Crank-Nicholson method in two respects. Firstly, there is a 
At 2 term multiplying the time difference term. This difference is unimportant because it does not alter 
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the accuracy of the scheme. Secondly, the cross-derivative terms are treated explicitly. If we ignore 
the cross-derivative terms, eq. (76) is as accurate as the Crank-Nicholson scheme, but is much easier to 
realize because it is easy to solve eq. (76) for / +1 . The explicit treatment of the cross-derivative terms 
lowers the accuracy and the stability, putting a limit on the maximum At that can be used. On the other 
hand, the implicit treatment of the other terms means that this method is far superior to the fully explicit 
method. 

We can compute / +1 from eq. (76) in a series of simple steps: 

4> k = h-Af k , 

e +i/2 = (i+^A v y\ k , 

* fc+i = V +i/2 , 

The inversion of the matrices is carried out using Gaussian elimination as described in appendix A. 

C. Example 

Let us consider a specific example relevant to lower hybrid current drive. The plasma consists of 
electrons and infinitely massive ions with Zi = 1. Electron-electron collisions are computed assuming a 
Maxwellian background using 

C mL e q- ( 41 >- Electron -ion collisions are given by eq. (37). The effect of 
the lower hybrid waves is modeled by a quasilinear diffusion coefficient given by eqs. (46) with D = 1, 
V\ = 3, and v 2 = 5. The electric field E is taken to be zero. Except for minor details this is the 
same example treated in the paper on lower hybrid current drive. 6 (The time normalization used in that 
paper differs from the one adopted here by a factor of two.) We take f(t = 0) = f m , v max = 10, 
M = N = 100, and At = 0.2. 

In studies of current drive, we are principally interested in the current density J, the rf power absorbed 
per unit volume by the plasma P, and their ratio J/ P. These are defined by eqs. (19) and (20) whose 
discrete forms read 

J= lnt(t;COSg) , (77) 
n 

^ M-l N 

P= -J2J2 2 nsm6 i+1/2 v*S wv , i+1/2 j AvA6, (78) 

where n is given by eq. (59). (These definitions include a 1/n factor, because the n is included in the 
normalizations for J and P.) 

The current is plotted as a function of time in fig. 4. With At = 0.5, the integration is unstable. The 
difference in the values of the current when the equations are integrated with At = 0.2 and At = 0.05 is 
about 0.1% of the final current. 

The steady-state solution for / is shown in fig. 5. This may be obtained by integrating the equation 
sufficiently long (until about t = 1000) with a fixed time step or else using the techniques described in 
sec. IX. (With this numerical method, the steady state is independent of At.) The plateau in the resonant 
region is clearly visible as well as the considerable perpendicular heating. Using eqs. (77) and (78), we 
have J = 5.754 x 10~ 2 , P = 4.011 x 10~ 3 , and J/P = 14.34. 

The flux plot for this case is given in fig. 6. This shows that the combination of rf diffusion and 
collisional scattering induces a perpendicular flux in the resonant region. Such flux plots are useful 
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in providing guidance for the analytic solution of this problem. 6 More extensive examination of this 
example can be found in the original paper 6 including projections onto the v\\ axis, slices at constant v±, 
etc. 

There are two possible sources of error in these results: errors arising from the finite boundary 
(i.e., because v max is finite) and errors arising from the finite mesh. The effect of the boundary can 
be determined by increasing v max to 20 (and increasing N to 200). In the steady state, this gives 
J = 5.759 x 10~ 2 , P = 4.012 x 10~ 3 , J/P = 14.35— changes of less than 0.1%. Thus for this 
particular problem, w max = 10 is adequate. 

The effect of the discrete spatial grid is found by varying Av and Ad. This we do by keeping 
"ma* = 10 varying M and N with M = N. Thus we have N = 10/ Av and AO = Awtt/10. The 
results for J and J/P are shown in fig. 7. We see that there is a lot of scatter in the data which arises 
because D w is discontinuous. As Av and A9 are varied, grid points (those on which the flux is defined) 
enter or leave the resonant region v\ < v\\ < v 2 . Each time this happens, there is a jump in J and P. 
As Av — > 0, J approaches its asymptotic value of about 5.6 x 10~ 2 and the convergence to this value 
is as Av. The finite mesh error in J with M = N = 100 is about 3%. This rate of convergence can be 
understood because J and P are exponentially dependent on v\ [J ~ exp(— \v\)] and v\ is determined 
only to within ±^Av. Thus the relative error in J and P is about exp(^viAv) - 1 rj ^viAv. This 
gives a relative error of 15% for v\ = 3, N = 100, v max = 10. The actual error is somewhat less than 
this because the boundary of the resonant region cuts across the grid lines and so v\ is in fact determined 
more accurately than was assumed here. Because J and P are both subject to the same error, the ratio 
J / P is more accurately given: convergence to the asymptotic value of 14.24 is as Av 2 and the value with 
M = N = 100 is in error by less than 1%. 

If instead we use the truncated electron-electron collision operator C t e / U e nc , the steady-state distribution 
function is rather similar to that shown in fig. 5. However, the flux plot fig. 8 shows a new eddy at low 
velocities due to the overall drift of the electrons with respect to the ions. (This plot is obtained with the 
same parameters as for fig. 6.) In this case, we find J = 7.092 x 10~ 2 , P = 4.294 x 10~ 3 , J/P = 16.52. 
The enhancement of the efficiency J/P comes about because momentum (and hence current) is no longer 
lost when tail electrons collide with bulk electrons. 

A check on the implementation of the C t e / U e nc is given by measuring the electrical conductivity. For 
Z l = 1, the exact conductivity is given by table I as J/E = 7.429 w 0.582 x 16^/2/tt. 21 Integrating 
the Fokker-Planck equation using the truncated collision operator with no rf Dq = and a small electric 
field E = 10~ 3 , the conductivity is J/E = 7.446, a 0.3% error. This small error is probably attributable 
partly to the finite mesh size (here we again took M = N = 100 and w max = 10) and partly to the 
finiteness of E (since there is a contribution to the current which varies as E 3 ). In contrast, if C^^., is 
used the conductivity is J/E = 3.772 a factor of two too small. 5 



IX. Steady-State Solution 
A. Statement of problem 

Often, we are only interested in the steady-state solution to the Fokker-Planck equation. Nearly 
always we must resort to an iterative method for obtaining the steady state. In that case we need some 
measure of how close we are to the steady state so that iteration may be stopped when this is small 
enough. The measure we shall employ is 



R = -4/int 
n V 



dt 



(79) 
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where the residue df/dt is given by eq. (60). Somewhat arbitrarily we use R = 10 9 as the convergence 
criterion. 

One obvious way of obtaining a steady state is to integrate the time-dependent solution as described 
in sec. VIII for a long time. This should be done with the largest time step consistent with stability. For 
the example shown in fig. 5, the convergence criterion is met at time t = 812. The largest time step 
that can be used is approximately 0.2; so that 4060 steps are required. The CPU time required to run the 
Fokker-Planck code on the Cray-1 is approximately 2 [is per mesh point per time step. Thus, achieving 
the steady state by this method takes about 80 s. This is rather expensive and it is therefore desirable to 
find faster methods. 

However, this method is quite effective when A x = 0. Then the numerical scheme is stable even 
if At is large. For example, for the electric field example discussed in sec. VIII in which D a = and 
E = 10~ 3 , we can take At = 1, and the convergence criterion is met after 220 steps. Here the integral 
portion of C t e / U ' nc , which is represented by the term h in eq. (71), is evaluated every tenth time step. 
The numerical method is stable for larger values of At. But, because the integration is less accurate, 
more steps are required to meet the convergence criterion. With large At the numerical solution tends to 
oscillate about the steady state. 

B. Chebyshev acceleration 

A significant improvement can be achieved by using a varying time step. Hewett et al. 32 describe an 
adaptive time selection for the alternating direction implicit method which speeds the convergence by a 
factor of two to three. Here we describe Chebyshev acceleration 31 which is a nonadaptive method for 
selecting varying time steps. We choose the time step Atk = tfc+i — tk according to 



where a, (3, and K are constants with a < (3 and K = integer. The advantage of this method is that by 
changing a few lines of code it can easily be incorporated into the alternating-direction-implicit method 
described in sec. VIII. A fixed time step is recovered in the special case a = (3 = 1/ At. 

Let us discuss the choice of the parameters in eq. (80). With K large, eq. (80) gives a series of K 
time steps (repeated periodically) varying from 1/a down to 1/(3. In the examples we consider, we take 
K = 20. Then the maximum time step is somewhat less than 1/a while the minimum time step is very 
close to 1/(3. In order to realize performance gains with this method we wish to pick the minimum time 
step comfortably within the stability threshold for the fixed-time-step method, while the maximum time 
step is considerably greater than the stability threshold. 

The method works because the long wavelength eigenmodes of the linear operator decay slowly but 
are stable with large At; on the other hand, the short wavelength modes decay rapidly but are only stable 
if At is small. Consider a particular cycle of K steps. During the initial large time steps, the long 
wavelength modes are efficiently damped (because At is large), but the short wavelength modes grow. 
This is followed by successively shorter time steps which damp the short wavelength modes. 

For the example shown in fig. 5, the stability threshold for At lies between 0.2 and 0.5. Thus we 
choose 1/(3 = 0.05 and 1/a = 1000. With K = 20 this gives a maximum time step of 3 1 .4, a minimum 
step of 0.05, and an average time step of 1.95. Since the average time step is about 10 times the largest 
time step that can be used in the fixed time step scheme, we expect convergence to be 10 times faster. 
Indeed this is the case. The convergence criterion is met after 400 steps at t = 790. This takes about 
8 s of CPU time. The variation of R with time is shown in fig. 9. This shows the growth of R during 
the large time steps followed by a drop in R as the instabilities are quenched during the small time steps. 



At k = 



2 



(80) 
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The overall decay of R with t closely matches that seen with a fixed time step. (This is contrary to the 
experience of Hewett et al. with their adaptive code in which the rates of decay are very different. 32 ) 

C. Runaway problem 

If the electric field is sufficiently large to produce runaways, i.e., E > u max , then as t — ► oo a 
steady state is reached which decays at the runaway rate 7 (assuming that a linear collision operator is 
employed). Because / and all its moments decay at the same rate, 7 is given from eq. (61) as 

mt 1 

which we will take to be the definition of 7 for all t. Thus we write 

/(v, t) = /'(v, t) cxp (- J 7 (0 dt'^j , (82) 

where 7 is given by eq. (81) and f'(t — > 00) is independent of t. If eq. (82) is substituted into eq. (71), 
we obtain 

df 

+ 7)/' = 0, (83) 

where for simplicity we set the inhomogeneous term h to zero. Because 7 is expressed as an integral 
over / eq. (81), it varies slowly and need not be evaluated very often. Thus eq. (83) may be regarded as 
a linear equation and solved in precisely the same way as eq. (71) (with h = 0) except that 7 must be 
subtracted from Vi+i/2.j+i/2 e Q- (73b). 

As an example, fig. 10 shows the steady-state distribution obtained by this method with Zi = 1, 
E = 0.06, M — N = 100, u max = 10, and electron-electron collisions given by C^ x . Since there 
is no rf diffusion term, there are no cross-derivative terms and the steady state is most easily obtained 
by taking a constant time step of At = 1. The runaway rate 7 is recomputed every ten time steps 
and the convergence condition R = 10~ 9 is met after 820 time steps. In the steady state, we have 
7 = 5.211 x 10~ 5 and J = 0.3133. These are close to the results obtained by Kulsrud et al, 5 namely 
7 = 5.411 x 10~ 5 and J = 0.3143. 

Again, it is important to explore the possible errors in these figures. Extending the boundary to 
u mM = 20 and doubling N to 200 gives 7 = 5.210 x 10~ 5 and J = 0.4514. While there is practically 
no change in 7, J is about 50% larger. This discrepancy arises because there is a large contribution to 
the total current by the runaways in the region 10 < v < 20. We can verify this by estimating the total 
current for an arbitrary w max on the basis of the results from v max = 10. For simplicity, assume that all the 
runaways are concentrated near v± = 0. From small 7 and in the limit t — > 00, the runaway distribution 
is independent of «||, so that /(^y » v t ) ~ ( / y/E)5(v±). The current obtained by integrating uy/ out to 
v = Wmax is then 

J(v max ) w Jbulk + \("f/E)v 2 max , 

where, using the data from v max = 10, we have Jbulk = 0.270. We can interpret Jbulk as the current 
carried by the bulk electrons and the other term as the current carried by the runaways. This now gives 
J(vma X = 20) = 0.444 which is within 2% of the observed value. The lesson from this exercise is that 
it makes little sense to quote the result for J when the runaway rate is appreciable because it depends 
strongly on v max . It is preferable to determine the bulk current since this is then weakly dependent on 
v max and has a physical interpretation. We have seen that v max — 10 is sufficiently large to give 7 and 
Jbulk accurately. 

In order to determine the effect of the finite mesh on the runaway results, we vary M and N with 
M = N and v max = 10. The results for 7 and J are shown in fig. 11. The asymptotic values are 
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7 = 5.185 x 1CT 5 and J = 0.31334. The errors in the values for M = N = 100 are 0.5% and 0.02%, 
respectively. The errors are considerably less than with the rf problem in fig. 7 and the convergence is 
much more regular (as Av 2 ). 

A disadvantage of solving for the decaying steady state of the distribution, eq. (83), is that S is 
no longer divergence free. This means that the stream lines cannot be plotted as contours of a stream 
function A, eq. (65). This can be remedied by injecting electrons at the origin to match the runaway loss 
of particles. Although this is a rather artificial problem, there is little error in the runaway rate provided 
that the runaway rate itself is small. We implement this procedure as follows: The loss of particles at 

V = Umax is 

n7 = flux 1 . 



We match this loss by a uniform radial flux at the origin 

v Q b v , i+ i/2,o - 2jtM 



which is chosen to give 



M-l 

i=0 



2ir sm0 i+1/2 v o l S Vyl+1/ 2,o A0 = n-f. 



(The product v 2 S V}i+ i/ 2 fi is finite even though £,,,1+1/2,0 is infinite.) From eq. (60), we see that this in- 
troduces a source term v 2 S v s+1/2,0/ { v \/ 2 i nto tne expressions for <9/,+i/2,i/2/c^. This is included 
as part of the inhomogeneous term h in eq. (71). The expressions for the stream function eqs. (65) require 
a slight modification to give 

2nv 2 

Aij = -7 H 2^ sm9 i , + i /2 S v ^ + i/2,j AO, 

n i'=0 
27rsin0, 3 ~ ] 



^ n sin vi y—- \ 

= Afi 2^ v j'+T-/^ s e,i,j'+i/2 Av, 



n 

j'=o 

where the integration constant has been chosen to given Aqj = —7 and Amj = 0. 

The flux plot computed by this method for the case shown in fig. 10, i.e., for Zi = 1, E = 0.06, 
M = N = 100, u max = 10, is shown in fig. 12. When computed in this way, the runaway rate is slightly 
lower 7 = 5.148 x 10~ 5 because a typical runaway particle has to be accelerated from v = instead of 
v = 1. The current J = 0.3127 is also lower. 



D. Other methods 

An infinite time step can be used if the Crank-Nicholson scheme, eq. (75), is modified so that f k+1 is 
used in place of \ (f k+1 + .f k )- Then, the steady state can be achieved in a single time step. Of course, this 
entails inverting the large matrix A (which is why we advocated using the alternating-direction-implicit 
method in preference to the Crank-Nicholson method). However, routines are available to perform such 
an inversion and they have been employed by O'Brien et aJ.. 33 An important feature of this method 
is the use of disk files to hold intermediate results. (Typically, the full matrix cannot fit into memory.) 
They report a CPU time of 35 s to invert the matrix arising from the discretization of the Fokker-Planck 
equation on a 300 x 100 grid with this time scaling as MN x min(M, TV). This method is therefore 
comparable (as far as CPU time goes) to the Chebyshev acceleration method. There are two potential 
drawbacks of this scheme: Firstly, there is a significant cost in I/O time with this method because of the 
use of disk files for storage. Secondly, the advantage of the method is reduced if the steady state cannot 
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be reached in a single time step. This is the case with the more complicated collision operators, because 
the matrix A is a function of time. 

Various iterative methods are available for obtaining a steady-state solution. 31 These are basically 
approximate methods of inverting the matrix A. Notable is Gauss-Seidel relaxation in which the elements 
of / are successively updated to achieve df/dt — at the point in question. In line relaxation, a 
whole line of elements (for example, j — const) is updated simultaneously (requiring the solution of 
a tridiagonal system of equations). Line relaxation gives the same convergence rate as Gauss-Seidel 
relaxation and may be vectorized if the even-numbered rows (j = even) are updated in one sweep 
followed by the update on the odd-numbered rows. 

The odd-even line relaxation method is extended with the successive-over-relaxation method where 
the over-relaxation parameter uj determines how much overshoot there is beyond the value of / which 
gives df/dt — 0. Unfortunately, these methods give results which are roughly the same as using fixed 
time steps. For the example shown in fig. 5, with the over-relaxation parameter set to u = 1.4, the 
convergence criterion is met after 5980 steps. (Compare this to the 4060 steps required in the fixed-time- 
step method. However, one relaxation step tends to be computationally less expensive than one step of 
the alternating-direction-implicit method.) For this example, the method becomes unstable with uj > 1.5. 

Although by themselves relaxation methods are not very useful for this problem, they are an important 
ingredient in the multigrid method. 34,35 In this method, the problem is solved at several different grid 
spacings (usually differing from each other by a factor of two). A few relaxation sweeps are carried 
out on the finest grid. Because relaxation is a local method, this is very effective at damping the short 
wavelength modes (with wavelength comparable to grid spacing). If relaxation is continued on the finest 
grid, convergence would become slower because longer wavelength modes would dominate the residue. 
However, in the multigrid method, the residue is transferred onto the next coarsest grid where relaxation 
methods are again efficient. This process continues recursively up to very coarse grids where either 
relaxation methods or direct solution methods can be used. 

This method has not been implemented for the Fokker-Planck equation. However, we can estimate 
the time required to obtain a steady state. Each relaxation step on the finest grid gives a reduction in R 
by about a factor of two. (The total work at the coarser grids is at most a multiple of the work on the 
finest grid.) In contrast, the mean reduction in R with the Chebyshev method is by 4% per step (see 
fig. 9). Thus the multigrid method will require about log(0.5)/ log(0.96) w 16 times fewer steps — an 
order-of-magnitude improvement over the Chebyshev method. 

X. Relativistic Treatment 

A. The Fokker-Planck equation 

Fokker-Planck methods have been used to study current drive by lower hybrid waves. In a fusion 
plasma, these waves will interact with electrons that travel at close to the speed of light. In such cases, it 
is necessary to reformulate the equation to include relativistic effects. The first change is that the electron 
distribution function is expressed in momentum rather than velocity space so that eq. (1) becomes 



where now the V = dj dp operator operates in momentum space, S w is the rf-induced flux in momentum 
space, and f e is normalized so that 




(84) 



S 
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In spherical coordinates we have 

Id Id 

V • S = -^—p 2 S p H —T^smOSe, 

p z op p sin 00 

where cos 6* = p\\/p. 

In addition, the forms of the collision term and the quasilinear diffusion term are altered. 

B. The relativistic collision operator 

The relativistic collision operator is given by Beliaev and Budker. 36 It can again be written as the 
divergence of a flux C(f a , fb) = — V • S" , where now we have 

s ./6 = ||lnA^| U(u) • (/a(p)^^ - A(p')^) d'p'- (85) 

The expression for U is rather complicated. 36 However, if either the test or the background species is 
weakly relativistic (p <C m a c or p' -C m^c), then U may be approximated by its nonrelativistic form 



U ( u ) = 3 , U = V a 



"6i 



where v s = p/m s j s is the velocity of species s, j s = (1 + p 2 /m^c 2 ) 1 / 2 is the relativistic correction 
factor, and m s is the rest mass. 

Despite the resemblance of eq. (85) to eq. (21), this collision operator cannot be readily expressed in 
terms of Rosenbluth potentials. However, considerable progress can still be made by working directly 
with eq. (85). We restrict our attention to electron-ion and electron-electron collisions. 

For collisions off infinitely massive ions, we can take the ions to be stationary v \ — > and evaluate 
the integrals to give 

C e /*(/ e (p)) =r e / e -^-l-|-sin0|-/ e (v), (86) 
v ' 2v e p sin oO oO 

where 

/b _ n b g 2 g 2 lnA°/ fc 
4ne 2 

(this differs by a factor of to 2 from the definition given in sec. II). 

For electron-electron collisions we start with the case of an isotropic background C(/ e (p), fe°\p))- 
The fluxes for this term are 19 



~,e/e _ n e/e 1 9f< 

> P de 

where 

47rr e/e / fP „,/2 



S% e = ~D%l--£, (87b) 



a) 



a -ne/e / r-p 12 roo -i 

*''•> =^(f«"«¥* ,+ f«"'<4 (88b) 

Fe d e = -^tr{j o P7 e (o) (po e Jr d P '+ j p p>f?\p>)2 Ve /c 2 dp>y m 
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These should be compared with their nonrelativistic counterparts eqs. (26) and (27). 
In the relativistic limit, the Maxwellian distribution eq. (18) becomes 37 

£ 



fem{p) 



n e 



47TTO2 C T e if 2 (e- 1 ) 



exp 



(89) 



where 

is the total electron energy, 



m e c 7 e 



9 = T e /m e c 2 = Te/511 keV, 

and K n is the nth-order modified Bessel function of the second kind. If we substitute (p) — f em (p) 
into eqs. (88), we obtain Fcl e /Dtpp — —v e /T e . Thus we find that f em annihilates the electron- 
electron collision term C(f em , f em ) = 0. The integrals in eq. (88) cannot be performed analytically 
with fe°^ (p) = fem (p) and so in the numerical code these are performed numerically. 
For the Maxwellian distribution eq. (89), we define a thermal momentum 

p te = V m eTe, 

a thermal velocity 



and a thermal collision frequency 



Vte = 



m, 



T e/e 



Pte 
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Forp » pte, the indefinite limits in the integrals in Eq. (88) can be replaced by oo, giving 

(90a) 
(90b) 
(90c) 

J-eVe 

These should be compared with eqs. (28). 

For a background which consists of just the first Legendre harmonic, the collision term is C(/em(p), 
fe 1 ^ (p) cos 6) . This is given by 19 

C(fe m (p)Je 1} (p) cos 9) _ 47rr e / e 



D e/e 
^ cpp 


pe/e V te 

3 ' 
V% 


n e/e 


= r e / e — ( i 

2v e \ 


pe/e 
cp 


-pe/e ^te 
TeVl 



femip) COS 9 

J m e / e (1) (p) 

I 7e 



lie 



+ 



P ' 2 fP(p')- 



7e V e 
P 2 7e 3 



_^ ( 4 7 - + 6)-i(4 7 --9yj 



+ 



i r°° 

■ I, 



P 2 Ye 3 

i 2 V e ( m e v'i 3 1 



dp' 



r(47 e 2 + 6)-i(4 7e 3 -9 7e ) 



-I- _i — _ 

p' 2 7 3 



3(47e 2 + 6) 



dp'}. 



(91) 



34 



[Compare with eq. (34).] The general solution of the linearized electron-electron collision operator 

C(f e J em ) + C(f em J e ) = is 

f e = (a + b-p + d£)f em , 

where a, b, and d are arbitrary constants. With a = d = and b = p |, this provides a useful check on 
Eqs. (88) and (91) and their computational realizations. 

In the example we show below, we use the electron-ion collision operator given by eq. (86) and the 
relativistic generalization of the truncated collision operator eq. (38) 

C t e /u„c(/e(p)) = C(/ e ( P ), femip)) + C(/ em (p), (p) COS 9) , (92) 

where the first term is given by eqs. (87) and (88) and the second term by eq. (91). 

C. Wave-particle interaction 

We saw in sec. V that the quasilinear diffusion operator had two principal ingredients: the wave- 
particle resonance condition, and the diffusion paths. Both of these are modified by relativistic effects. 
The wave-particle resonance condition becomes 

w - k\\v e \\ - nO e /7 e = 0, 

where Cl e = q e B/m e is the rest-mass cyclotron frequency. Translating this into momentum space gives 

tjJ\J\ +p 2 /mlc 2 — — nO e = 0. 

This modification of the resonance condition is important in the consideration of current drive by electron 
cyclotron waves. 39 

The diffusion paths are again given by surfaces of constant energy in the wave frame. The expression 
for the energy in a frame moving at (w/fc||)p|| is 

= s - (oVfcj|)pj| 

The diffusion paths are, therefore, given by 

£ — {w/k\\)p\\ = const. 

These paths are parallel to the vector 

^-« e ||)p±+«e±P||- 

This should be compared with the vector a„ defined in sec. V. The paths are ellipses or hyperbolae in 
momentum space depending on whether co/kn is less than or greater than c. 25 

For lower hybrid waves, we have n — and the diffusion is in the parallel direction. We, therefore, 
generalize eqs. (46) by incorporating the modified resonance condition to read 

D w = Au(rl,P||)P||P||, (93a) 

where 
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D. Example 

To illustrate the relativistic effects we show in fig. 13 the steady-state distribution function obtained 
by integrating the Fokker-Planck equation with electron-electron collisions given by C^ nc eq. (92) and 
electron-ion collisions given by C e l % eq. (86) with Z{ = 1. The quasilinear diffusion term is given by 
eq. (93) with D = l,vi = 0.4c, and v 2 = 0.7c. (Except for the perpendicular profile of D w , this is the 
same as the example given in ref. 19.) The integration is carried out with M = N = 100 and p max = 20. 
We normalize all momenta to pt e , velocities to ptejm e (not Vt e ), the current density to n e q e pte/iTi e , the 
power density to n e p1 e v te /m e , etc. Again we are principally interested in the current and the power 
dissipated. These are defined by 



where int is the generalization of eq. (57) to momentum space, n = int 1. [Compare these expressions 
with eqs. (77) and (78).] In the steady state, we find J = 3.732 x 10~ 3 , P = 1.256 x 10~ 4 , and 
J IP = 29.72. 

Again a useful benchmark is provided by the electrical conductivity. In the limit E — > this is 
correctly given if C t e / U e nc is employed. With E = 10" 3 , Z, = 1, 6 = 0.01, M = N = 100, and 
Pmax = 10, we find J/E = 7.307, which differs from the true value of 7.291 by about 0.2%. Values of 
the conductivity for various values of and are tabulated in table II. 

XI. Adjoint Method 

A. Introduction and example 

We have considered here techniques for solving the Fokker-Planck equation with an added quasi- 
linear diffusion term. This tends to be an expensive operation because the addition of the quasilinear 
diffusion term greatly increases the parameter space to be scanned. For example, the study of lower hy- 
brid current drive 6 included the results of some 50 runs with different values of v\ and u 2 . Even so, no 
systematic study was made of the dependence on the parameters D n and Zj. 

However, the amount of work can be drastically reduced using the adjoint method. This was intro- 
duced by Hirshman 40 for the study of beam-driven currents. Later, Antonsen and Chu 41 used it to study 
rf-driven currents. 

To illustrate the method, we will outline the analysis given by Antonsen and Chu. 41 The method 
begins by assuming that f e is close to a Maxwellian f em so that the linearized electron-electron collision 
operator C^ e eq. (31) can be used. The quasilinear diffusion term is taken as a given. As pointed out in 
sec. V, the Fokker-Planck equation then becomes an inhomogeneous equation, whose linear operator is 
independent of the wave drive. Two further assumptions are made, namely that E = and that a steady 
state has been reached. (Neither of these assumptions is necessary and they have been relaxed in ref. 42.) 
The Fokker-Planck equation is then 



int ( t; e cos 9) 




n 

M-l N 



i=0 j=0 



C(/ e (v))^C 1 < :/ e (/ e (v))+C e / l (/ e (v)) 




(94) 
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where we have inserted the Chapman-Enskog-Braginskii energy loss term to ensure that eq. (94) has 
a solution (i.e., to ensure that the Fokker-Planck equation reaches a steady state). Taking the energy 
moment of this equation, and noting that the collision operator is energy conserving, we find the equation 
fori)I. ill 

where P is given by eq. (20). 

The straightforward approach is now to solve eq. (94) for a particular S w , determine the electron 
distribution f e , and hence find the rf-driven current. The adjoint method gives a way of computing the 
current without having to find f e . Consider first the "adjoint" problem 

C(/em(«)x(v)) = -q e v\\f em (v), (95) 

where we require that f em x have zero density and zero energy. This is the Spitzer-Harm equation for the 
perturbed distribution in the presence of an electric field E = T e v\\. Let us multiply eq. (95) by f e /f em 
and integrate over velocity. This gives 



J = ~ J(fe/fe m )C(f emX ) d 3 V, 



where J is the current carried by the electron distribution f e . Now we utilize the self-adjointness of the 
linearized collision operator 

J i>C(f emX ) d 3 V = J xC(fem^) d 3 V, 

together with eq. (94) for C(/ e ) to give 

J= f S„(v).V X (v)d 3 v. (96) 



Equation (96) is the desired expression for the current. The quantity \ serves as the Green's function for 
the current J. The current drive efficiency is given by 

P J m e S w ■ vd 3 v 

B. Solving the adjoint equation 

In order to apply this method, we must determine \ by solving eq. (95). Because x(v) consists of only 
the first Legendre harmonic x' 1 ^ (v) cos 9, this equation reduces to a one-dimensional integro-differential 
equation, 

where D e J v % and D e J g e are given by eqs. (30), and l e l e is defined by 

= C(f em (v), f em (v)x {1) (v) cos 9) 

(X ( ^ fem(v) COS 8 

[seeeq. (34)]. 
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In general, eq. (98) must be solved numerically. This is done by constructing the partial differential 
equation by setting the left-hand side of eq. (98) equal to dx^/dt. The resulting equation is integrated 
in time with arbitrary initial conditions until a steady state is reached. The integration is carried out 
in the domain < v < v max and the boundary conditions are taken to be x ( v = 0) = and 

d 2 x {1) (v = v max )/dv 2 = o. 

Approach to the steady state is accelerated by treating the first three terms in eq. (98) fully implicitly; 
i.e., in order to compute [x^'ty + At) — x^'(t)]/ At, we evaluate these terms at t + At. This means that 
very large time steps can be used. The integral term I e / e {x^) is treated explicitly and is reevaluated 
at each time step. The resulting difference equations form a tridiagonal matrix which can be solved by 
Gaussian elimination. 

Because the adjoint equation is the same as the equation for the perturbed distribution in the presence 
of a weak electric field, we can solve eq. (98) to obtain values of the electrical conductivity which is 
defined by 

i = S ^? /em(wMv)d3v= ^ / v3 f^)x (1) (v)dv. 

This procedure was carried out using the method outlined above with w max = 15v te , Av — 0.001ut e , and 
At = lOOO/t'te. Because we are only working with a one-dimensional equation, it is possible to use a 
much finer mesh than with two-dimensional problems and so obtain results which are effectively "exact." 
The results are summarized in table I where we have also included the results from use of approximate 
collision operators. The same technique is easily generalized to relativistic plasmas using the collision 
operator given in sec. X. This gives the relativistic corrections to the conductivity which are given in 
table II. 

When the adjoint method is applied to more complicated situations (e.g., including a dc electric field), 
a two-dimensional equation must be solved. We can then use many of the techniques for the solution of 
the Fokker-Planck equation, which have been presented in the preceding sections. 

C. Discussion 

Let us assess the work involved in utilizing the adjoint method. Once the adjoint equation has been 
solved, the current and the efficiency are immediately given in terms of S w by eqs. (96) and (97). Instead 
of having to solve the Fokker-Planck equation afresh for every form of S w , a couple of velocity integrals 
over S„, suffice to give the important quantities. The parameter space that must be scanned in order to 
give a complete understanding of the physics is greatly reduced. The adjoint method does not give the 
electron distribution f e nor the rf-induced flux S^,. On the other hand, a crude estimate of S„, gives 
an accurate estimate of the efficiency because eq. (97) involves the ratio of two integrals over S„,. An 
effective way to use this method within a ray-tracing code would be to determine S w from a solution of 
the one-dimensional Fokker-Planck equation 1 and to use this to determine both J and P from eqs. (96) 
and (20). The code thereby benefits from an accurate determination of the current drive efficiency while 
the high computational costs of integrating the two-dimensional Fokker-Planck equation are avoided. 

Because the current drive efficiency is determined by a single function x, it is possible to ask questions 
not readily answerable from numerical solutions of the Fokker-Planck equation. Examples are: What is 
the asymptotic form for the efficiency as the wave phase velocity becomes large? What is the maximum 
possible efficiency for a particular class of waves? 

Besides determining the current, the adjoint method can be adapted to give other moments of the 
electron distribution by changing the right-hand side of eq. (XI). This can then give, for example, the 
perpendicular energy of the electrons, bremsstrahlung radiation, etc. This method has been used to 
determine the current-drive efficiency in a relativistic plasma. 19 Recent developments of the method 42 
allow the determination of arbitrary moments of f e (not just the current J), and the determination of the 
time development of such moments. These have been applied to the study of rf current ramp-up. 43 
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XII. Conclusions 



In the last fifteen years, Fokker-Planck codes have gone from esoteric programs developed by a few 
researchers which could only be run on a few machines to widely available tools used by a large number 
of physicists on many different computers. This has been due to the large increase in computer power 
available to the average physicist and the pioneering efforts of Killeen et aJ.. 2,3 

In this paper, I have given a detailed description of a particular implementation of a code to solve the 
Fokker-Planck equation with emphasis on a particular application, namely current drive by lower hybrid 
waves. There are many other implementations of this code that have been applied to a large variety of 
interesting problems. My goal has been to illustrate the main numerical problems by means of concrete 
examples. The methods presented here cover the major numerical problems that are encountered in all 
Fokker-Planck codes. 

There are two areas which still require attention. Firstly, improved methods for obtaining the steady- 
state solution of the Fokker-Planck equation are needed. Here the multigrid method offers the best 
promise for substantial savings over the other methods described in this paper. Secondly, the adjoint 
methods outlined in sec. XI should be extended and applied to a wider range of problems. Ray-tracing 
codes still need to be modified to accept the results of these calculations. 
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Appendix A. Numerical Techniques 

In this appendix, various fragments of code are shown. A two-dimensional Fokker-Planck code is 
ideally suited to a vector processing machine like the Cray-1. However, care must be taken to order the 
loops correctly, otherwise they will not vectorize. 

The first example is the computation of the current int(wcos6 | ) eq. (57). This illustrates the rather 
peculiar way in which FORTRAN code must be written in order to take advantage of the Cray-l's archi- 
tecture. 44 We assume that the arrays and variables given in table III have been initialized as indicated. 

dimension temp(0 : iy — 1) 
do 1 i = 0, iy — 1 
temp(i) = 0.0 

1 continue 

do 3 j = 0, jx - 1 
do 2 i — 0, iy - 1 
temp(i) = temp(i) + x(j)**3 * f(h.j) 

2 continue 

3 continue 

do 4 i = 0, iy — 1 
temp(i) — sn(i) * cn(i) * temp(i) 

4 continue 
cur = 0.0 

do 5 i = 0, iy — 1 
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cur = cur + temp(i) 
5 continue 

cur = 2.0 *pi* dx * dy * cur 

The important point is that the inner loop (with label 2) vectorizes. This would not happen if the order 
of the loops were reversed. There is no particular advantage in taking the computation of x(j)**3 out 
of the inner loop since the CFT compiler does this automatically. The only loop that the compiler treats 
inefficiently is the last one. In fact, we replace this by a call to the OMNILIB routine ssum. 

The second example is computing the integral part of the truncated collision operator C(/ m (t>), 
(v) cos 9) I cos 9 eq. (34). Here again it is easy to arrange so that most of the code vectorizes. 4 The 
computation of this term is then relatively inexpensive compared with the other computations. 

dimension s0{— 1 : jx — 1), s3(0 : jx), s5(0 : jx),fl (0 : jx — 1) 
do 1 j = 0,jx — 1 

/JO') = o.o 

1 continue 

do 3 i = 0, iy - 1 
do 2 j = 0, jx - 1 
f 1 CO = f 1 U) + 1-5 * dy * sn(i) * cn(i) * f(i,j) 

2 continue 

3 continue 

do 4 j = 0,jx — 1 
s0(j-l) = dx*fl(j) 
s3(j + 1) = s0(j - 1) * x(j)**3 
s5(j + 1) = s3(j + 1) * x(j)**2 

4 continue 

s3(0) = 0.5* s3{l) 
s5(0) = 0.5* s5(l) 
do 5 j = 1, jx — 1 

s3(j) = sS {j - 1) + 0.5 * (s3(j) + s3(j + 1)) 

s5 (j) = s5(j - 1) + 0.5 * (s5(j) + a5 (j + 1)) 

5 continue 

s0(jx- 1) = 0.5* s0{jx- 2) 
do 6 j = jx-2,0, -1 
s0(j) = s0(j + 1) + 0.5 * (s0(j) + s0(j - 1)) 

6 continue 

do 7 j = 0, jx — 1 
cl(j) = (s5(j)/5.0-s3(j)/3.0)/x(j)**2 
+ + s0(j) * (x(j)**2/5.0 - 1.0/3.0) * x(j) 

cl (j) = 4 * pi * fm(j) * (/I (j) + cl (j)) 

7 continue 

All the loops vectorize with the exception of the indefinite integration loops (with labels 5 and 6). Most 
of the time is spent in the inner loop 2 during the computation of /W eq. (15). 

Finally, we consider vectorized Gaussian elimination. This subroutine performs Gaussian elimination 
for the tridiagonal system of equations 

to give Xij for < i < n, < j < m. The coefficients satisfy a^j = c„_ij = 0. A substantial fraction 
of the running time of the Fokker-Planck code is spent in this subroutine. When implemented for a single 



40 



system of equations m = 1, this leads to "vector dependencies" which inhibit vectorization. The solution 
is to solve the m systems in parallel with j being the index for the inner loops. In the subroutine below, 
it is assumed that all the matrices are the same size, that the spacing in memory between Xij and 
(the solution direction) is ns, and that the spacing between Xij and Xij+i (the vectorizing direction) is 
ms. This subroutine uses x and y as temporary storage; thus the initial data in y are destroyed. 

subroutine solve(a;, ns, n, ms, m, a, b, c, y, dt) 
dimension ir(0 : ms — 1, : m — 1), y(Q : ms — 1, : m — 1), 
+ a(0 : ms — 1, : m — 1), 6(0 : ms — 1, : m — 1), 

+ c(0 : ms - 1,0 : m - 1) 

dt§ = 0.5 * rft 
do ^ i = 0, n - 1 
ia = ns * (i — 1) 
ife = ns * i 
do 1 j — 0, m — 1 
rfen = 1.0/(1.0+ cfcg * (6(i6, j) + a(ib,j) * y{ia,j))) 
x(ib,j) = {y{ib,j) — dt^ * a(ib,j) * x(ia,j)) * den 
y(ib,j) = —dt2 * c(ib,j) * den 

1 continue 

2 continue 

do 4i = n- 2,0,-1 
i6 = ns * i 
ic = ns * (i + 1) 
do 5 j = 0, m - 1 
x{ib,j) = y(ib,j) *x(ic,j) +x(ib,j) 

3 continue 

4 continue 
return 
end 

There are a couple of tricky points here. Firstly, we use nonstandard indexing into the arrays. The element 
Xij is accessed by the array element x(ns * i, j). If ms = 1, then ns * i will generally exceed the upper 
bound m.s — 1 on the first dimension of the arrays. This type of array indexing may cause problems with 
compilers that perform bounds checking. Secondly, we have utilized the fact that a(0) = c(n — 1) = 
and assumed that an arbitrary (possibly undefined) number multiplied by zero will give zero. If this is 
not the case, the i = and i = n — 1 iterations in the loop with label 2 will have to be split off from the 
rest of the loop and treated separately. 

This subroutine is sufficiently general to be used for both the matrix inversions required in imple- 
menting eq. (76). Assuming that all the matrices are dimensioned by, for example, 

dimension /(0 : iyl — 1, : jxl — 1) 

then the inversions are obtained by 

call solve(xia, iyl,jx, 1, iy, ax, bx, cx,phi, dt) 
call so\ve(xib, 1, iy, iyl,jx, ay, by, cy, xia, dt) 
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Tables 



TABLE I. The electrical conductivity for various values of the ion charge Zi and for various electron- 
electron collision operators. The conductivities are normalized to n e ql/m e v te . 



Collision operator 


Zi = 1 


Zi = 2 


Zi = 5 


Zi = 10 


linearized 


7.429 


4.377 


2.078 


1.133 


drifting 


6.331 


3.876 


1.932 


1.084 


Maxwellian 


3.773 


2.824 


1.660 


0.998 


high-velocity 


2.837 


2.310 


1.489 


0.938 



TABLE II. The electrical conductivity of a relativistic plasma for various values of the ion charge Zi 
and for various electron temperatures. The conductivities are normalized to n e q 2 / m e v te and the electron 
temperatures are given in terms of 9 = T e /m e c 2 . 



e 


Zi = 1 


Zi = 2 


Zi = 5 


Zi = 10 


0.0 


7.429 


4.377 


2.078 


1.133 


0.01 


7.291 


4.275 


2.019 


1.097 


0.02 


7.160 


4.180 


1.963 


1.064 


0.05 


6.807 


3.928 


1.821 


0.979 


0.1 


6.317 


3.590 


1.636 


0.872 


0.2 


5.575 


3.102 


1.383 


0.729 



TABLE III. Meaning of FORTRAN variables and arrays. 



FORTRAN name 


meaning 


dx 


Av 


dy 


A6 


dt 


At 


jx 


N 


iy 


M 


xg(j) 




x{i) 


V 3 + l/2 


yg(i) 


0i 


y(i) 




cg{i) 


cos 9i 


cn(i) 


COS i+ i/2 


sg(i) 


sin#i 


sn(i) 


sin^+i/2 


pi 


7T 


f(hj) 


/i+l/2,j+l/2 


fm{j) 


/m,j+l/2 


cur 


int(u cos 9) 


fl(j) 


/ (1) (^ + l/ 2 ) 


cl (j) * cn(i) 


C(f m (v),fW(v)co S 6)\. +1/w/2 


ax(i,j) 


a vd+l/2,j + l/2 


ay(i,j) 


a 0,i+l/2,j + l/2 


phi(i,j) 


4+1/2J + 1/2 
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v ll 

FIG. 1. The cylindrical and spherical coordinate systems. 




FIG. 2. The relation between the resonance condition for quasilinear diffusion u> — k\\V\\ — nVL e = and 
the diffusion path (v — (cj/fc||)v||) 2 = const. 
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t 



FIG. 4. The current as a function of time for Zi = 1, f(t = 0) = f m , and rf diffusion given by eqs. (46) 
with D = 1, vi = 3, and v 2 = 5. Here we have M = N = 100, At = 0.2, and u max = 10. 
Electron-electron collisions are computed using CJ^. 
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FIG. 5. The steady-state distribution for the case shown in fig. 4. The contour levels are / = (2-7r)~ 3 / 2 x 
cxp[— i(j'/5) 2 ] for j = integer. This gives equally spaced contours for a Maxwellian distribution with 
spacing 6v = |. The resonant region is shown. 
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FIG. 7. The current J (a) and the efficiency J/P (b) as functions of Av. The parameters are the same as 
for fig. 4 except that M and N are allowed to vary with M = N. The plots show the results from runs 
with N varying between 100 and 350 in steps of 5 and between 350 and 500 in steps of 50. 
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FIG. 8. The flux plot when C t e / U e nc is used. The parameters are otherwise the same as for fig. 6. 



R 




t 



FIG. 9. R as a function of time when Chebyshev acceleration is applied to the example shown in fig. 
Here 1/(3 = 0.05, l/a = 1000, K = 20, and f(t = 0) = f m . The convergence criterion R < 10~ 9 
met after 400 steps at t = 790. 
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FIG. 10. The steady-state distribution in the presence of a dc electric field. Here we have Zi = 1, 
E = 0.06, M = N = 100, w max = 10 and electron-electron collisions are computed using C^ax- The 
contour levels are the same as for fig. 5. 
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FIG. 11. The runaway rate 7 (a) and the current J (b) as functions of Av. The parameters are the same 
as for fig. 10 except that M and TV are allowed to vary with M = N. The plots show the results from 
runs with N varying between 50 and 300 in steps of 10 and between 300 and 500 in steps of 50. 
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FIG. 12. The flux plot for the runaway problem. This illustrates the same case as shown in fig. 10 except 
that a source of particles is introduced at the origin to balance the runaway loss 7 = 5.148 x 1CT 5 . One 
set of contour levels is 0.1j(j + |) for j = integer and —10 < j < 10 (these give the stream lines that 
run away and the outermost stream lines that encircle the central eddy). The other set of contour levels is 
2 x 10~ 4 (j + |) for j = integer and j > (these are the innermost stream lines about the eddy). 




FIG. 13. The steady-state distribution for Z, = 1, 6 = 0.01 (T e = 5.11 kcV), and rf diffusion 
given by eq. (93) with Do = 1, v\ = 0.4c, and v 2 = 0.7c. Here we have M — N = 100 and 
Pmax = 20. Electron-electron collisions are computed using C t e / U e nc . The contour levels are chosen 
to be / = y/Qexp[—^yi + Q(j /3) 2 /Qi\/[A-n:K2{Q~ 1 )] for j = integer which give equally spaced 
contours for a relativistic Maxwellian with spacing Sp = i. [For O = 0.01 we have ^(O 1 ) = 
1.019-^/71-/2-^/0 cxp(— 1/6).] The resonant region is shown. 
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